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Analytical  solutions  for  diffuse  interface  propagation  are  found  for  two  recently  developed 
Landau  potentials  that  account  for  the  phenomenology  of  stress-induced  martensitic  phase 
transformations.  The  solutions  include  the  interface  profile  and  velocity  as  a  function  of 
temperature  and  stress  tensor.  An  instability  in  the  interface  propagation  near  lattice  insta¬ 
bility  conditions  is  studied  numerically.  The  effect  of  material  inertia  is  approximately 
included.  Two  methods  for  introducing  an  athermal  interface  friction  in  phase  field  models 
are  discussed.  In  the  first  method  an  analytic  expression  defines  the  location  of  the  diffuse 
interface,  and  the  rate  of  change  of  the  order  parameters  is  required  to  vanish  if  the  driving 
force  is  below  a  threshold.  As  an  alternative  and  more  physical  approach,  we  demonstrate 
that  the  introduction  of  spatially  oscillatory  stress  fields  due  to  crystal  defects  and  the  Pei- 
erls  barrier,  or  to  a  jump  in  chemical  energy,  reproduces  the  effect  of  an  athermal  thresh¬ 
old.  Finite  element  simulations  of  microstructure  evolution  with  and  without  an  athermal 
threshold  are  performed.  In  the  presence  of  spatially  oscillatory  fields  the  evolution  self¬ 
arrests  in  realistic  stationary  microstructures,  thus  the  system  does  not  converge  to  an 
unphysical  single-phase  final  state,  and  rate-independent  temperature-  and  stress-induced 
phase  transformation  hysteresis  are  exhibited. 

©  2009  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Phase  field  or  Ginzburg-Landau  (GL)  models  are  widely  used  for  the  simulation  of  various  first-order  solid-solid  phase 
transformations;  see  books  Salje  (1991)  and  Toledano  and  Toledano  (1987).  In  this  paper  we  will  focus  on  martensitic  or  dif¬ 
fusionless  transformations  (Ahluwalia  et  al.,  2003;  Artemev  et  al.,  2001;  Curnoe  and  Jacobs,  2001a, b;  Jacobs  et  al.,  2003;  Jin 
et  al.,  2001;  Levitas  and  Preston,  2002a, b;  Levitas  et  al.,  2003;  Levitas  and  Lee,  2007;  Lookman  et  al.,  2003a, b;  Rasmussen 
et  al.,  2001;  Seol  et  al.,  2003;  Shenoy  et  al„  1999;  Wang  et  al.,  2001).  Deformation  of  the  crystal  lattice  of  the  austenite, 
A,  the  high-temperature  phase,  into  the  martensite,  M,  the  low  temperature  phase,  can  be  described  by  the  transformation 
strain  tensor  st  (also  called  the  Bain  strain  or  spontaneous  strain).  The  relative  symmetries  of  the  A  and  M  crystal  lattices 
implies  the  existence  of  a  finite  number  n  of  crystallographically  equivalent  variants  of  martensite.  All  martensitic  variants 

Mj,  i=l,2, _ n,  have  the  same  components  of  the  transformation  strain  tensor  in  their  respective  crystallographic  bases.  A 

list  of  components  of  transformation  strain  tensors  for  transformations  between  various  crystal  lattices  are  given,  for  exam¬ 
ple,  in  Bhattacharya  (2004)  and  Pitted  and  Zanzotto  (2002).  The  phase  field  approach  describes  both  stress-  and  tempera¬ 
ture-induced  phase  transformations  in  a  unified  framework. 


*  Corresponding  author.  Tel.:  +1  515  294  9691:  fax:  +1  801  788  0026. 
E-mail  address:  vlevitas@iastate.edu  (V.l.  Levitas). 

0749-641 9/$  -  see  front  matter  ©  2009  Elsevier  Ltd.  All  rights  reserved. 
doi:l  0.1 01 6/j.ijplas.2  009.08. 003 


396 


V.I.  Levitas  et  al.  / International  Journal  of  Plasticity  26  (2010)  395-422 


In  the  phase  field  framework  the  evolution  of  a  multi-connected  martensitic  microstructure  is  described  by  a  thermody¬ 
namically  consistent  set  of  kinetic  equations,  the  Ginzburg-Landau  equations,  for  the  order  parameters.  Order  parameters 
are  similar  to  internal  variables  in  continuum  thermodynamics  (Valanis,  1996).  However,  in  the  phase  field  approach  these 
internal  variables  describe  material  instabilities,  such  as  the  instabilities  of  a  crystal  lattice  responsible  for  solid-solid  phase 
transformations,  twin  and  dislocation  nucleation,  melting,  fracture  and  so  on.  Some  theories  of  martensitic  phase  transfor¬ 
mations  employ  order  parameters  related  to  transformation  strains  (Artemev  et  al.,  2001 ;  Levitas  and  Preston,  2002a, b;  Lev¬ 
itas  et  al.,  2003;  Levitas  and  Lee,  2007;  Seol  et  al.,  2003;  Shenoy  et  al„  1999;  Wang  and  Khachaturyan,  1997;  Wang  et  al„ 
2001),  while  the  order  parameters  in  other  models  are  the  components  of  the  strain  tensor  responsible  for  lattice  instability 
(Curnoe  and  Jacobs,  2001a,b;  Jacobs  et  al.,  2003;  Lookman  et  al.,  2003a, b;  Rasmussen  et  al.,  2001).  The  thermodynamic  (Lan¬ 
dau)  potential  is  typically  a  polynomial  in  the  order  parameters  with  multiple  minima  (as  in  Fig.  2)  corresponding  to  the 
various  phases.  The  phase  with  the  deepest  minimum  is  the  stable  phase,  while  other  minima  correspond  to  metastable 
phases;  all  minima  are  separated  by  potential  barriers.  The  thermodynamic  potential  also  includes  gradient  terms.  Solutions 
of  the  time-dependent  Ginzburg-Landau  equation,  which  describes  the  evolution  of  the  order  parameters,  are  generally 
comprised  of  regions  corresponding  to  local  minima  of  the  potential  (stable  or  metastable  phases)  separated  by  diffuse  sta¬ 
tionary  or  moving  interfaces  where  the  gradient  energy  is  localized;  it  represents  the  interface  energy.  The  key  advantage  of 
the  phase  field  approach  is  that  the  computation  of  the  microstructure  evolution  proceeds  without  the  additional  effort  re¬ 
quired  to  track  multiple  interfaces.  However,  the  standard  phase  field  method  does  not  encode  the  microphysics  governing 
interface  propagation. 

As  we  demonstrated  in  Levitas  and  Preston  (2002a, b),  all  previous  Landau  potentials  did  not  account  for  the  typical  fea¬ 
tures  of  stress-strain  curves  for  martensitic  phase  transformations,  e.g.,  in  shape  memory  alloys.  In  our  papers  (Levitas  and 
Preston,  2002a, b;  Levitas  et  al.,  2003)  we  developed  several  polynomial  Gibbs  (Landau)  potentials  for  the  description  of 
multivariant  stress-  and  temperature-induced  martensitic  phase  transformations  under  general  three-dimensional  loading. 
These  potentials  were  designed  by  requiring  that  they  describe  the  experimentally  observed  features  of  martensitic  phase 
transformation  in  shape  memory  alloys  and  steels,  specifically,  a  constant  or  weakly  temperature  dependent  transformation 
strain  tensor  and  stress  hysteresis,  and  transformation  at  non-zero  tangent  elastic  moduli.  They  include  all  temperature- 
dependent  thermomechanical  properties  of  the  austenite  and  martensitic  variants  and  describe  phase  transformations  be¬ 
tween  austenite  and  martensitic  variants  and  between  martensitic  variants  for  arbitrary  crystal  structures.  These  potentials 
are  based  on  order  parameters  related  to  transformation  strain  rather  than  total  strain.  We  do  not  know  how  to  derive  a  sim¬ 
ilar  potential  in  terms  of  order  parameters  related  to  the  total  strain. 

In  Levitas  et  al.  (2003),  analytic  solutions  of  the  one-dimensional  time-independent  Ginzburg-Landau  equations  for  our  2- 
3-4  and  2-4-6  polynomial  potentials  were  obtained.  Solutions  include  martensitic  (M)  and  austenitic  (A)  critical  nuclei,  and 
diffuse  M-A  and  M-M  interfaces.  The  widths  and  energies  of  the  nuclei  and  interfaces  were  found  as  functions  of  the  ther¬ 
modynamic  driving  force,  the  gradient  energy  coefficient,  and  a  parameter  that  characterizes  the  stability  of  A.  Static  micro¬ 
structures  in  a  finite  sample,  their  stabilities,  and  physical  interpretations  were  studied  in  Levitas  et  al.  (2006a).  Combined 
surface  and  size  effects  and  a  barrierless  nucleation  mechanism  were  analyzed  in  Levitas  et  al.  (2006b).  Dynamic  problems 
were  treated  in  Idesman  et  al.  (2008). 

In  this  paper,  interface  propagation  kinetics  is  incorporated  in  the  models  developed  in  Levitas  and  Preston  (2002a, b)  and 
Levitas  et  al.  (2003).  We  consider  both  one-dimensional  (ID)  and  two-dimensional  (2D)  cases.  In  the  ID  case,  we  obtain  and 
analyze  both  analytical  and  numerical  solutions  of  the  time-dependent  Ginzburg-Landau  equations  for  both  A-M  and  M-M 
interface  propagation.  Both  2-3-4  and  2-4-6  polynomial  potentials  are  used.  For  the  A-M  interface  an  exact  solution  for  the 
interface  profile  and  velocity  is  obtained  for  negligible  inertia  (mass  density),  in  which  case  the  stress  is  the  same  in  both  A 
and  M,  i.e.,  aA  =  <rM  =  a.  Analytical  relationships  between  the  interface  velocity,  the  driving  force  for  the  phase  transforma¬ 
tion,  and  a  parameter  that  characterizes  the  lattice  stability  of  A  are  obtained  and  analyzed.  For  non-zero  mass  density, 
cta^o-m.  To  a  good  approximation  the  driving  force  for  interface  propagation  depends  on  aA  and  <rM  only  through  their  aver¬ 
age.  Thus,  inertial  effects  can  be  approximately  taken  into  account  by  replacing  the  homogeneous  stress  a  in  the  zero-inertia 
A-M  interface  solution  by  a  =  ( aA  +  <xM)/2. 

The  propagation  of  M-M  and  A-M  interfaces  is  studied  numerically.  For  the  M-M  interface,  we  simulate  the  case  where 
the  temperature  equals  the  phase  equilibrium  temperature  and  an  austenitic  region  appears  between  the  martensitic  vari¬ 
ants.  We  present  several  numerical  solutions  illustrating  instabilities  in  M-M  and  A-M  interface  propagation.  When  stresses 
reach  and  exceed  the  values  corresponding  to  lattice  instability,  homogeneous  or  surface-induced  nucleation  may  occur  in 
addition  to  interface  propagation,  and  the  nucleation  interacts  with  the  interface  propagation. 

In  2D,  the  coupled  Ginzburg-Landau  and  quasi-static  equations  of  linear  elasticity  theory  are  solved  using  the  finite  ele¬ 
ment  method  (FEM).  The  coupled  evolution  of  microstructures  and  stress  fields  in  square  samples  is  studied  numerically. 

Despite  significant  success  in  modeling  microstructure  formation  in  Artemev  et  al.  (2001),  Curnoe  and  Jacobs  (2001  a, b), 
Jacobs  et  al.  (2003),  Jin  et  al.  (2001 ),  Levitas  and  Preston  (2002a, b),  Levitas  et  al.  (2003),  Levitas  and  Lee  (2007),  Lookman  et  al. 
(2003a, b),  Rasmussen  et  al.  (2001),  Seol  et  al.  (2003),  Shenoy  et  al.  (1999),  Wang  and  Khachaturyan  (1997),  and  Wang  et  al. 
(2001),  and  here,  the  phase  field  approach  has  a  major  drawback:  it  does  not  include  an  athermal  resistance  to  interface  mo¬ 
tion.  This  resistance  is  analogous  to  dry  friction  in  classical  mechanics.  Because  of  this  athermal  resistance,  interface  prop¬ 
agation  occurs  only  if  the  driving  force  for  the  corresponding  phase  transformation  exceeds  a  rate-independent  threshold  K. 
The  athermal  resistance  is  responsible  for  hysteresis  in  the  temperature  or  the  rate-independent  part  of  the  stress,  and  en¬ 
ergy  dissipation. 
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In  our  numerous  simulations  (some  of  them  are  presented  below)  for  a  single  crystal  with  homogeneous  stresses  at  the 
boundary  we  found  that  complex  martensitic  microstructures  appear  that  are  similar  to  those  observed  experimentally; 
however,  they  eventually  evolve  into  a  single  phase.  Similar  results  have  already  been  reported  in  the  literature  (see,  for 
example,  Jacobs  et  al.,  2003).  Consequently,  microstructural  formation  in  samples  with  stress-free  surfaces  (temperature-in¬ 
duced  phase  transformations),  with  homogeneous  stresses  at  the  boundaries  (as  in  the  experiments  in  Abeyaratne  et  al. 
(1996)),  or  with  zero  stresses  at  selected  surfaces  (for  example,  uniaxial  tension-compression  or  torsion),  cannot  be  mod¬ 
eled.  Even  for  periodic  boundary  conditions,  the  final  microstructure  is  sometimes  a  single  variant  (Jin  et  al„  2001 ;  Kerr 
et  al.,  1999).  In  contrast,  kinematic  constraints,  e.g.,  due  to  polycrystallinity  or  prescribed  displacements  at  the  boundary, 
promote  stationary  multivariant  microstructures  (Jacobs  et  al.,  2003;  Jin  et  al„  2001;  Rasmussen  et  al.,  2001).  However, 
an  athermal  threshold  exists  in  any  case,  and  if  taken  into  account  it  would  change  the  microstructure  evolution,  path  depen¬ 
dence,  and  energetics,  especially  under  cyclic  loading.  While  the  necessity  of  introducing  an  athermal  threshold  in  phase 
field  modeling  has  been  recognized  for  a  long  time,  we  are  not  aware  of  any  successful  attempts  to  do  so. 

An  athermal  threshold  is  included  in  all  mesoscale  models  (Auricchio  et  al.,  2007;  Boyd  and  Lagoudas,  1996;  Ghosh  and 
Olson,  1994;  Grujicic  et  al.,  1985;  Levitas,  1994, 1995, 1998,  2000a, b;  Levitas  et  al.,  1999,  2002c;  Levitas  and  Ozsoy,  2009a, b; 
Lim  and  McDowell,  2002;  Pan  et  al.,  2007;  Peng  et  al„  2008;  Thamburaja  and  Anand,  2002)  that  include  kinetic  equations  for 
sharp  interfaces  or  for  product  phase  concentration.  It  is  also  included  in  our  recent  mesoscale  phase  field  model  (Idesman 
et  al.,  2005;  Levitas  et  al.,  2004)  since  neglect  of  the  gradient  energy  term  (which  is  reasonable  at  the  mesoscale)  results  in  a 
threshold  value  for  the  driving  force  for  interface  motion.  However,  the  introduction  of  an  athermal  threshold  in  the  tradi¬ 
tional  nanoscale  phase  field  approach  (where  gradient  energy  cannot  be  neglected)  is  not  straightforward  because  there  is  no 
equation  for  the  interface,  only  evolution  equations  for  the  order  parameters.  Thus,  Vedantam  (2006)  introduced  a  sophis¬ 
ticated  kinetic  coefficient  in  the  Ginzburg-Landau  equation  that  is  singular  for  zero  rate  of  change  of  the  order  parameter  17, 
which  allows  one  to  obtain  fj  =  0  for  non-zero  driving  force.  However,  the  introduction  of  such  a  threshold  in  the  kinetic 
equation  for  the  order  parameters  arrests  certain  unphysical  intermediate  configurations  (e.g.,  particular  critical  nuclei); 
therefore  the  system  does  not  converge  to  a  realistic  microstructure  consisting  of  austenite  A  and  martensitic  variants  M, 
divided  by  moving  or  fixed  diffuse  interfaces.  Note  that  a  similar  problem  exists  in  the  phase  field  theory  of  dislocations 
(Hu  et  al.,  2004;  Wang  et  al.,  2001 ).  In  Wang  et  al.  (2001 ),  a  periodic  thermodynamic  potential  in  terms  of  order  parameters 
was  introduced  for  dislocations.  The  potential,  which  was  called  a  ‘Peierls  potential’,  was  supposed  to  represent  the  Peierls 
barrier  to  interface  propagation  due  to  discreteness  of  the  crystal  lattice,  but  it  was  shown  in  Hu  et  al.  (2004)  using  an  ana¬ 
lytical  solution  for  the  dislocation  that  there  was  in  fact  no  athermal  threshold.  Similarly,  the  analytical  expressions  for  inter¬ 
face  velocity  in  our  original  GL  model  (Levitas  and  Preston,  2002a, b;  Levitas  et  al.,  2003)  do  not  include  an  athermal 
threshold.  More  generally,  there  is  no  threshold  in  any  GL  model  because  the  gradient  energy  term  renders  it  non-local,  thus 
interface  points  are  coupled  through  the  order  parameter  and  interface  motion  occurs  at  any  non-zero  driving  force.  Since 
there  is  no  essential  difference  between  GL  theories  for  phase  transformations  or  dislocations,  a  prescription  for  introducing 
an  athermal  threshold  for  interface  motion  applies  equally  well  to  dislocations. 

An  athermal  threshold  appears  in  phase  field  theory  when  the  latent  heat  of  transformation  is  taken  into  account  (Ngan 
and  Truskinovsky,  1999).  Allowance  for  microinertia  related  to  the  strain  tensor  also  produces  a  threshold  (Theil  and  Levitas, 
2000).  A  discrete  model  can  exhibit  athermal  hysteresis;  see  Kressea  and  Truskinovsky  (2003).  The  quasi-continuum  approx¬ 
imation  of  a  discrete  model  can  also  exhibit  some  hysteresis,  but  it  is  much  smaller  than  for  the  initial  discrete  model  (Kres¬ 
sea  and  Truskinovsky,  2003). 

None  of  the  above  approaches  solves  the  problem  of  athermal  friction  for  interfaces  or  dislocations.  Indeed,  for  slow  inter¬ 
face  motion,  the  hysteresis  due  to  latent  heat  release  or  internal  inertia  is  too  small  to  account  for  measured  values  of  K.  In 
fact,  apart  from  the  Peierls  barrier,  the  origin  of  the  athermal  threshold  K  is  the  interaction  of  the  moving  interface  or  dis¬ 
locations  with  the  long-range  stress  fields  of  point  and  line  defects,  and  various  boundaries  (e.g.,  twin  and  tilt  boundaries) 
(Ghosh  and  Olson,  1994;  Grujicic  et  al.,  1985;  Kocks  et  al.,  1975).  The  most  compelling  evidence  to  support  this  claim  is  the 
observed  proportionality  of  the  athermal  threshold  and  the  rate-independent  hysteresis  to  the  yield  strength,  which  char¬ 
acterizes  the  types,  densities,  and  distributions  of  the  point,  line,  and  boundary  defects  that  limit  dislocation  motion;  see 
Levitas  (1997,  1998,  2004)  for  high  pressure  phase  transformations,  and  Ghosh  and  Olson  (1994)  and  Levitas  et  al. 
(2002c)  for  martensitic  transformations  in  steel.  Also,  it  is  known  that  the  phase  transformation  hysteresis  is  changed  along 
with  the  defect  microstructure  by  thermomechanical  treatment  (Hornbogen,  1999). 

In  this  paper,  we  propose  two  different  schemes  for  introducing  an  athermal  interface  friction  in  phase  field  models,  and  we 
show  that  each  provides  a  realistic  description  of  interface  propagation.  In  the  first  method,  an  analytic  expression  defines  the 
location  of  the  diffuse  interface,  and  the  rate  of  change  of  the  order  parameters  in  a  neighborhood  of  the  interface  is  required  to 
vanish  if  the  driving  force  is  below  some  threshold.  In  that  case,  the  interface  is  subject  to  an  athermal  friction  similar  to  that 
experienced  by  sharp  interfaces  in  micromechanical  models,  and  the  value  of  K  is  the  only  required  information.  This  method 
works  well  in  the  1 D  case  but  should  be  checked  for  2D  and  3D  problems.  As  an  alternative  approach,  we  demonstrate  that  the 
introduction  in  our  phase  field  model  of  spatially  oscillatory  stress  fields  (due  to  the  Peierls  barrier  and  various  defects),  or  of  a 
jump  in  chemical  energy,  AG",  reproduces  the  effect  of  an  athermal  threshold  on  interface  propagation.  In  the  presence  of  spa¬ 
tially  oscillatory  fields,  experimentally  observed  microstructures  self-arrest  before  the  system  can  converge  to  a  single-phase 
final  state,  and  rate-independent  temperature-  and  stress-induced  phase  transformation  hysteresis  are  exhibited.  These  fields 
do  not  affect  the  evolution  of  an  intermediate  microstructure,  such  as  a  critical  nucleus,  to  an  A— M,  or  M,— Mj  microstructure 
divided  by  diffuse  interfaces.  Also,  some  experimentally  observed  microstructures  appear  in  phase  field  simulations  with 
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oscillatory  fields  that  do  not  follow  from  GL  energy  minimization.  A  similar  approach  can  be  applied  to  other  phase  transfor¬ 
mations  -  reconstructive,  ferroelectric,  ferroelastic,  and  magnetoelastic  -  as  well  as  to  dislocation  motion.  The  incorporation  of 
spatially  oscillatory  fields  in  phase  field  models  dramatically  improves  the  fidelity  of  numerical  simulations  of  the  develop¬ 
ment  of  martensitic  microstructures.  However,  a  new  problem  arises:  the  necessity  of  finding  realistic  oscillatory  fields  cor¬ 
responding  to  a  given  defect  structure;  this  is  just  as  important  but  even  more  challenging  than  finding  a  proper  Landau 
potential.  We  also  discuss  the  introduction  of  an  athermal  threshold  in  a  model  for  thermally  activated  dislocation  motion. 
The  results  of  many  FEM  simulations  of  microstructure  evolution  with  an  athermal  barrier  are  presented. 

Note  that,  practically  speaking,  the  phase  field  approach  is  limited  to  nanoscale  samples  because  the  width  of  the  diffuse 
interface  is  of  the  order  of  magnitude  of  1  nm  and  it  has  to  be  resolved  with  several  finite  elements.  In  the  given  paper,  all 
material  parameters  for  the  thermodynamic  potential  were  determined  from  the  results  of  molecular  dynamic  simulations 
(see  Levitas  and  Preston,  2002b).  At  such  small  scales  (for  example,  in  nanofilms),  one  usually  neglects  nucleating  crystal 
defects.  Then,  there  are  no  stress  concentrators  to  promote  nucleation,  and  nucleation  occurs  close  to  the  lattice  instability 
at  stresses  of  the  order  of  magnitude  of  10  GPa.  These  values  correspond  to  molecular  dynamic  simulations  while  stresses 
measured  in  macroscopic  experiments  are  below  1  GPa.  This  situation  is  similar  to  the  case  of  fracture  and/or  plastic  flow: 
stresses  that  cause  fracture  and  plastic  flow  in  defect-free  nanoscale  volumes  correspond  to  the  theoretical  strength,  which  is 
of  the  order  of  magnitude  of  10  GPA;  the  engineering  strength  and  yield  strength  for  materials  with  defects  are  1-2  orders  of 
magnitude  lower.  In  contrast,  interface  propagation  can  occur  close  to  the  thermodynamic  equilibrium  conditions,  and  we 
operate  with  stresses  of  several  100  MPa. 

Some  authors  (e.g.,  Vedantam,  2006)  calibrate  the  parameters  in  the  thermodynamic  potential  using  experimental  data 
for  macroscopic  samples.  With  these  parameter  values,  interface  widths  and  energies  are  calculated  to  have  values  much 
larger  than  those  observed.  The  problem  is  that  the  parameter  values  determined  from  macroscopic  data  incorporate  the 
effects  of  mesoscale  defects  and  should  not  be  used  at  small  length  scales.  The  microscale  phase  field  approach  is  based 
on  essentially  different  concepts  (see  Levitas  et  al.,  2004;  Idesman  et  al.,  2005). 

2.  Interface  propagation  in  one  dimension 

In  this  section  we  provide  the  background  and  establish  notation  for  our  analytical  and  numerical  studies  of  A-M  and  M- 
M  interface  propagation  in  one  dimension.  All  transformations  can  be  described  with  a  single  order  parameter.  The  order 
parameter,  ij,  is  a  function  of  the  coordinate  x  in  the  direction  of  propagation  n;  n  n  =  1.  Numerical  simulations  are  carried 
out  in  a  rectangular  parallelepiped  in  an  arbitrary  three-dimensional  homogeneous  stress  field  a\  the  corresponding  normal 
and  shear  stresses  on  the  surface  are  denoted  cr,-  (t  =  1 , 2)  and  t  (Fig.  1 ).  The  transformation  strain  which  transforms  the  crys¬ 
tal  lattice  of  A  into  the  lattice  of  M  is  the  invariant  plane  strain  Et  =  ly(mn  +  nm )sign(ij)  +  enn,  where  y  is  the  shear  strain  in 
direction  m  in  the  habit  plane  with  normal  n  (direction  of  interface  propagation)  and  normal  strain  e  along  n;  the  faces  of  the 
parallelepiped  are  orthogonal  and  parallel  to  n  and  m.  The  transformation  strain  as  a  function  of  order  parameter  is  de¬ 
scribed  by  the  following  equations: 

Ete  =  Et(p6(r])  and  sf4  =  (1) 
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Fig.  1.  Schematic  of  solution  of  the  Landau-Ginzburg  equation  in  ID.  The  curve  labled  r\(x)  represents  the  solution.  The  crystal  lattice  transforms  from 
austenite  to  martensite  by  an  invariant  plane  strain. 
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where 


cp6{t])  =  aif/2  + (3  -  a)i]4  +  (a- 4)r]6/2  and  <pA{ij)  =  arj2  +  (4  -  2a)r]3  +  (a  -  3)rj4, 


(2) 


and  a  is  a  material  parameter,  0  ^  a  ^  6.  Here  and  later  the  subscripts  4  and  6  refer  to  the  2-3-4  and  2-4-6  potentials, 
respectively.  The  2-4-6  and  2-3-4  potentials  derived  in  Levitas  et  al.  (2003)  are 

C6  =  Si,72[1  -  (4  -  P)if/2  +  (3  -  P)t]4/ 3]/2,  (3) 

C4=s1f72[l  -(6-P)i]/3  +  (4-P)i]2/4},  (4) 

S,  :=  A  -  aa  :  £r,  S2  :=  12(AG"  -  a  :  £t),  P:=S2/S,.  (5) 


The  order  parameter  satisfies  -1  C  t]  1  for  the  2-4-6  potential  and  0  <  r]  ^  1  for  the  2-3-4  potential.  AG"  is  the  difference 
between  the  thermal  parts  of  the  Gibbs  energies  of  M  and  A.  Note  that  elastic  strains  are  neglected  in  Eqs.  (3)-(5),  but  can  be 
trivially  included.  The  minima  of  the  potential  G4  are  at  rj  =  0  (A),rj  =  1  (M),  and  for  G6  the  minima  are  at 
i]  =  0  (A),  i]  =  1  (M+),  i]  =  -1  (M_),  independent  of  a  and  temperature  9  (Fig.  2).  The  quantity  s,  characterizes  the  stability 
of  the  austenite  lattice;  Si  =  0  corresponds  to  the  loss  of  A  lattice  stability,  i.e.,  to  the  disappearance  of  the  A  minimum.  For 
both  potentials,  the  thermodynamic  driving  force  for  the  M  — >  A  transition,  i.e.,  G(l) — G(0),  is  equal  to  s2/ 12.  However,  the 
rate  of  the  M  — ►  A  transformation  is  controlled  not  just  by  the  driving  force  but  also  by  the  height  of  the  potential  barrier 
separating  the  A  and  M  minima.  The  location  and  height  of  this  barrier  are  simple  rational  functions  of  Si  and  s2 ;  in  the  limit 
s2  — >  0  of  small  driving  force  the  barrier  height  is  G4  =  (l/16)Si,  G6  =  (2/27)Si. 

Adding  a  gradient  energy  term  to  the  Landau  free  energy  yields  the  Ginzburg-Landau  energy  CGL  =  C  +  p(di]/dx)2  and 
Get  =  I,  GCLdx.  The  time-dependent  Ginzburg-Landau  (TDGL)  equation  follows  from  the  assumption  that  the  generalized  flux 
dij/dt  is  proportional  to  the  generalized  force  - SGa/Sr y. 


di]  _  ,  SGcl 
dt  ^  Si] 


(6) 


Here  SGCl/Si]  is  the  functional  derivative  of  Ga  with  respect  to  i]\  X  >  0  and  p  >  0  are  the  kinetic  and  gradient  energy  coef¬ 
ficients  with  dimensions  of  volume  I  energy-time  and  energy  I  length,  respectively.  The  GL  equation  for  the  2-3-4  potential  can 
be  written 


-nin- 1)  '7- 


di]  _  Err] 
dF~dx* 


T  =  (4si  -  s2)7t,  x!  = 


2s, 


4s,  -  s2  ’ 


4s,  -  s2 


2  P 


-x, 


and  the  GL  equation  for  the  2-4-6  potential  is 


dt]  _  Sri] 


—2  =  — 2  -  rjin2  -  T )  [if  - 

dt'  dx'2  n  1 


3s,  -  s2  ’ 


t'  =  (3s,  -  s2)Xt,  x'  = 


3s,  -  s2. 


2  P 


(7) 


(8) 


Note  that  Eqs.  (7)  and  (8)  are  nonlinear  diffusion  equations,  and  as  such  possess  solitonic  solutions  representing  propagating 
interfaces. 


Fig.  2.  Plot  of  Landau  potential  versus  order  parameter  for  the  potentials  G4  and  G6. 
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The  equations  of  motion  read 


da 


9 u  dr 


d2v 
'  dt2 


where  p  is  the  mass  density,  u  and  v  are  the  displacements  in  the  directions  n  and  m,  respectively,  a  =  n  a 
stress  (er  =  <7,  in  Fig.  1 )  and  t  =  n  a  m  is  the  shear  stress;  it  follows  that  a  s t  =  ere  +  vysign(t/). 


O) 

n  is  the  normal 


3.  Analytical  kink  solutions:  propagating  A-M  diffuse  interfaces 

3.1.  Zero  inertia 

3.1.1.  Interface  profile 

The  neglect  of  inertia  (mass  density)  in  the  equations  of  motion,  Eq.  (9)  means  that  the  stresses  a  and  t  are  constants; 
hence,  Sj  =  A  -  a[ae  + Tysign(r/)]  and  s2  =  AG"  -  as  -  vy sign(ri)  are  nearly  constant  provided  the  temperature  change 
accompanying  the  phase  transformation  is  small.  We  consider  the  case  where  the  phase  is  A  as  x  — ►  —  oo  and  M  as 
x  — >  +oo.  When  A  and  M  are  in  thermodynamic  equilibrium,  then,  G(-oo)  =  G(A)  =  G(+oo)  =  G(M)  =  0  and  s2  =  P  =  0.  The 
solutions  of  the  static  version  of  Eq.  (6)  read  (Levitas  et  al.,  2003) 

>lT(x)  =  [l  +  exp(-v/V^(x-x0))]  \  f/fcM(x)  =  [l  +  exp  (~y/2sA/P(x  -  x0))]  V2.  (10) 

The  transformation  strain  profiles  <p4[t/(x)]  and  tpe[r](x)]  and  the  ij(x)  profiles  are  shown  in  Fig.  13  in  Levitas  et  al.  (2003)  and 
schematically  in  Fig.  1.  They  smoothly  connect  the  austenitic  (t/  =  0)  and  martensitic  (tj  =  1)  regions. 

Propagating  interface  solutions  of  the  TDGL  equation,  Eq.  (6),  for  our  2-3-4  and  2-4-6  potentials  for  constant  stresses  a 
and  r  (zero  mass  density  in  the  equations  of  motion,  Eq.  (9))  can  be  obtained  by  generalizing  the  forms  of  and  r/^  given 
in  Eq.  (10).  In  the  case  of  the  2-3-4  potential  we  write 

»?4M(0  =  f1  +exP  +  ,  C  =  x  -  c4t.  (11) 


Here  c4  is  the  interface  velocity  and  F  is  a  function  to  be  determined.  Substituting  Eq.  (1 1 )  in  Eq.  (6)  and  requiring  the  coef¬ 
ficients  of  powers  of  exp(-  sjs\  /fit/)  and  exp(F(f))  to  vanish,  one  obtains  simultaneous  equations  for  dF/df  and  dF2 /dij2.  These 
equations  are  quadratic  in  dF/df;  signs  are  chosen  so  that  cs2  >  0  (i.e„  when  c  <  0  and  the  M  phase  grows  to  the  left,  then 
s2  <  0)  and  dF/df  =  0  for  s2  =  0.  The  F  derivatives  are  given  by 


dF/dC^/l-^1^2,  dF2/d£2  = 


4  \As] 


pa 


P 


2y/Jj  ’  4X 

Since  dF/d£  is  constant,  then  dF2/dt) 2  =  0,  which  determines  the  interface  velocity  c4.  Finally,  we  obtain 
1 


1+e  2Zf 


~(x-c4t) 


c4  =  .  =  2^/pots2sign(s2), 

V4si  -  s2 


(12) 


(13) 


where  a  :=  P/(4  -  P)  =  s2/(4sj  -  s2)  (a  =  0  for  thermodynamic  equilibrium,  a  =  1  when  M  loses  its  stability,  and  oc  =  -1 
when  A  loses  its  stability).  For  the  2-4-6  potential  we  have 


l2p  /S2 

Ce  V  3  \73s,  -  s2 


2  X^/2fJu.s2 

V3  s/3^a 


sign(s2). 


(14) 


3.1.2.  Interface  velocity 

In  contrast  to  the  sharp  interface  approximation  in  which  the  interface  velocity  depends  only  on  the  driving  force  for  the 
phase  transformation,  in  the  phase  field  approach  there  is  also  an  explicit  dependence  on  Si  (or  a)  which  characterizes  the 
relative  stabilities  of  the  phases.  However,  in  Eq.  (13)  for  c4,Si  appears  only  in  the  combination 

4s,  -  s2  =  4(A(0)  -  a<T :  st)  -  u(AGe  -  a  :  sX  (15) 

though  in  Eq.  (14)  for  c6.s,  appears  in  the  combination  3s,  -  s2  =  (3/4)(4s,  -  s2)  -  s2/4.  Using  the  approximate  relations 
(Levitas  and  Preston,  2002a, b) 

A  =  A0(8-9C),  AG6  =  Ao(0  -  0e)/3, 


(16) 


V.I.  Levitas  et  al.  / International  Journal  of  Plasticity  26  (2010)  395-422 


401 


where  0e  is  the  equilibrium  temperature  for  stress-free  A  and  M,  0C  is  the  temperature  of  A  loss  of  stability  for  zero  stress,  and 
Aq  —  -3A s  with  As  for  transformation  entropy,  we  find 

4s,  -s2  =  A0(6e  -  9C)  +  (3  -  a)a  :  st.  (17) 

Remarkably,  the  temperature  dependencies  of  A(0)  and  G°  have  cancelled,  and  since  a  is  approximately  equal  to  3  (Levitas 
and  Preston,  2002a, b),  we  set  a  =  3,  thereby  eliminating  the  dependence  on  stress  as  well;  therefore 

4s,-s2=A0(6e-Bc)  =  A  (18) 


is  a  constant.  With  these  approximations  it  follows  that  a  and  the  interface  velocities  in  both  potentials  are  functions  of  only 
the  driving  force,  as  in  the  sharp  interface  approximation 


a  = 


s2 

4A ' 


c4  = 


2 s2VP 

2  VT 


2Xs2\j2fi 

V3>\J\2A-s2 


(19) 


in  this  approximation,  a  is  proportional  to  s2  and  has  no  additional  dependency  on  stresses.  For  the  2-3-4  potential,  the 
interface  velocity  is  proportional  to  the  driving  force  for  the  phase  transformation,  but  for  the  2-4-6  potential  the  depen¬ 
dence  is  nonlinear;  it  is  stronger  than  linear  for  M  — >  A  phase  transformation  and  weaker  than  linear  for  A  — >  M  phase 
transformation. 


3.3.3.  Interface  width 

The  AM  interface  width  is  defined  in  Levitas  et  al.  (2003)  by 

a  am  fdtpMo  iv1 


V  dc  )r 


which  results  in 


AT  = 


2(21  -  5a  +  Y4 


32(a  -  6)3(1  la3  -  81(9  +  Y4)  -  5a2(24  +  Y4)  +  a{ 486  +  39Y4))  V  4s>  “  S2’ 

4 


(72  -  15a  +  Y6) 


3/3 


6  128\/2(a  -  6)3(-27a2  -  24(24  +  Y6)  +  5a(48  +  Ys))  V  3si  “  S2’ 

Y4  :=  \/81  -  30a  +  5a2,  Y6  :=  ^576-2400  +  3302. 

Both  A™  and  AgM  are  complicated  functions  of  a  but  they  are  accurately  approximated  by  the  polynomials 
1.88  ^  p6(a)  =  1.88  +  0.179a  +  0.00065a2  -0.0035a3  ^  2.386, 

4.822  <  p4(a)  =  -0.056(a  -  3)2  +  5.334  s;  5.334. 

Then  we  obtain 


A4m  =  p4(a) 


P 


'  4s,  -  s2 
in  the  approximation  of  Eq.  (19) 


=  P4(a)\  KM=Pe(a) 

V  S2 


3/3 


3s,  -  s2 


=  2ps(a) 


3/3a 


(3  -  a)s2' 


1271  -s2 


(20) 


(21) 


(22) 


(23) 


(24) 


The  interface  thickness  for  the  2-3-4  potential  is  independent  of  the  driving  force,  but  for  the  2-4-6  potential  it  increases 
(decreases)  with  driving  force  for  the  M^A(A->M)  phase  transformation. 


3.3.4.  Evaluation  of  A 

To  obtain  an  order  of  magnitude  estimate  of  the  kinetic  coefficient  A  we  take  a  =  3  and  s,  =  0  (loss  of  A  stability),  for 
which  s2  =  -4 A,  and  approximate  the  interface  velocity  by  the  shear  wave  velocity  cs.  From  Eq.  (19)  we  obtain 

cs  =  2l4\[fA  =  24  \J2pAf3.  (25) 

Using  the  following  data  for  NiAl  (Levitas  and  Preston,  2002a, b;  Levitas  et  al.,  2003;  Miracle,  1993) 

710  =  4.40  MPa  IC1,  /?  =  2.59  x  10~10  N.  0e  =  215K,  0C  =  -183K,  p  =  5850  kg/m3, 

G=71.5GPa,  and  cs  =  sjG/ p  =  3496  m/s, 


(26) 
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we  obtain 

/L,  =  \f2/3Xs  =  2596  m2/N  s.  (27) 

Note  that  p  =  2.59  x  1CT10  N  corresponds  to  the  equilibrium  width  of  1  nm  for  a  M-M  interface  (Levitas  and  Preston, 
2002a, b;  Levitas  et  al„  2003).  For  a  M-M  interface  width  of  0.3  nm  (which  corresponds  to  an  interatomic  distance)  the  gra¬ 
dient  energy  and  kinetic  coefficients  are  respectively  decreased  and  increased  to 

P  =  2.33  x  10~niV  and  24  =  fa/Vl  =  8653.5  m2/N  s.  (28) 

3.2.  Inertial  effects  on  interface  propagation 

The  results  of  the  previous  section  are  valid  in  the  limit  of  zero  mass  density,  which  implies  constant  stresses  er  and  t.  The 
actual  stress  variation  and  the  corresponding  solutions  for  p  are  governed  by  the  coupled  equations  of  elasticity  theory  and 
Eqs.  (6)  and  (9).  An  analytic  solution  is  impossible,  so  we  have  developed  an  approximation  scheme. 

We  consider  a  diffuse  planar  interface  separating  A  at  x  -oo  from  M  atx-i  +oo;  the  unit  normal  to  this  plane,  n,  is 
directed  from  A  to  M.  The  stress  tensor  a,  which  has  the  normal  component  a  and  in-plane  shear  component  t,  varies  from 
cta  at  x  — >  —  oo  to  the  prescribed  value  <rM  at  x  — »  oo. 

The  driving  force  for  M  — >  A  is  s2/12  =  AG"  -  a  :  et,  where  a  ■.  et  is  the  transformation  work  per  unit  volume.  Consider  an 
interface  separating  A  at  -oo  from  M  at  +oo  that  is  propagating  at  constant  speed  +c.  The  transformation  work  is  obtained  by 
integration 

/°°  dp 

|29) 

where  £  =  x  -  ct.  We  define  interface  profiles  fa  and  fe  for  the  stress  and  transformation  strain 

a{Q  =  a»  +  (a*-*»)fAQ, 

*(0  =  «£(£)• 

The  profiles  satisfy  fae  ~  1  for  f  — >  +oo  and/CT  E  ~  0  for  £  — >  -oo.  lffa  =f,  then  the  transformation  work  is  independent  of  the 
profile 


Wt  =  ^  (<rA  +  <rM)  :  £t,  (31) 

but  in  general  fa^ft-  In  this  case  we  approximate  the  profiles  by  quartic  polynomials  constrained  to  satisfy /^(O)  =  0  and 
fa.eiCo)  ~  1  (£0  is  the  interface  thickness),  and  //.{(0)  =/^e(f0)  =  0.  The  resulting  one-parameter  interface  profiles  are 

faAO  =  «*AM o)2  +  (4  -  2aff,e)(f/Co)3  +  ~  3)(C/£0)4  (32) 

which  is  just  <p4(£/£0)  with  a  replaced  by  a„( .  The  transformation  work  is 

Wt  =  +  aM)  :  £t  +  ^ (aff  -  a£)  ■  ^ (ffA  -  ffM)  :  £t.  (33) 


Note  that  Wt  is  independent  of  the  interface  thickness.  The  contribution  of  (1/2)(<ta  -  cM)  :  st  to  Wt,  which  arises  from  the 
difference  between  the  stress  and  strain  profiles,  is  suppressed  by  roughly  an  order  of  magnitude  (3/35)  relative  to 
(l/2)(c-A  +  o-M)  :  £t-  To  a  good  approximation,  Wt,  and  therefore  the  driving  force  depends  on  the  stresses  in  the  A  and  M  only 
through  their  average:  a  ■=  (o-A  +  o-M)/2;  hence,  Wf  =  a  ■.  e,.  Thus,  in  this  approximation,  inertial  (finite  mass  density)  effects 
are  taken  into  account  by  replacing  the  constant  stress  a  by  a.  This  is  consistent  with  the  sharp  interface  approach  (Abey- 
aratne  and  Knowles,  1993;  Kondaurov  and  Nikitin,  1986).  Accordingly,  we  utilize  Eqs.  (13)  and  (14)  with  a  — >  a  for  the  gen¬ 
eral  dynamic  case. 

Combining  the  Hadamard  compatibility  condition  [v]  =  -c[F]  ■  n  with  the  jump  condition  (o-M  -  <rA)  ■  n  =  — pc[v]  corre¬ 
sponding  to  the  equation  of  motion  (9),  one  obtains  for  sharp  interface 

(itm  -  <rA)  ■  n  =  pc2 [F]  ■  n.  (34) 

Here  [a]  =  aM  -  aA,  v  is  the  particle  velocity,  and  F  is  the  deformation  gradient.  Neglecting  the  jump  in  elastic  strain  across  the 
interface,  one  obtains  [F]  =  ymn  +  enn  for  invariant  plane  strain.  Then  Eq.  (34)  reads 

<jM  -  erA  =  pc2e,  tm  -  ta  =  pc2y.  (35) 

In  the  definitions  of  Si  and  s2,  a  :  e,  =  as  +  ty.  Considering  aM  and  tm  as  prescribed  values,  Eqs.  (35)  and  (13)2  (or  Eq.  (14)2) 
constitute  three  equations  in  the  three  unknowns  crA,  ta  (or  equivalently  a  and  f ),  and  c.  The  solutions  for  both  potentials  are 
given  by 
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Z4(jM  -  yeZ3  iM  -  fiZ7 
2(y2  +  £2)Z, 


Z5tm  -  y£Z3<7M  -  yZ7 
T=  2(72  +  62)2, 


(36) 


2  Z3(ecrM  +  7TM)+Z7 

p(y2  +  £2)Z,  '  1  J 

with 

Z,  =  m[l(y2  +  £2)22p,  Z2  =  n-a  +  Z,,  Z3  =  n-a  +  2Z,, 

Z4  =  (n-3)(272  +  £2)  +  272Z,,  Z5  =  (n  -  3)(y2  +  2e2)  +  2fi2Z,,  Z6  =  4  -  AG>  +  2Z,),  ,,Q, 

1/2 

Z7  =  Z6  -  (jZ6  +Z3(£< tm  +  7Tm)]2  -  4Z,Z2(£<tm  +  yrM  -  AC")2)  , 

where  (m,n)  =  (18,3)  for  the  2-3-4  potential  and  (m,n)  =  (16,4)  for  the  2-4-6  potential.  The  sign  of  c  is  chosen  from  the 
condition  cs2  >  0,  i.e  that  the  stable  phase  grows.  To  choose  the  physically  relevant  solution  when  solving  the  quadratic  we 
require  a  =  aM  and  i  =  rM  in  the  limit  of  zero  inertia  ( p  =  0).  The  above  equations  describe  the  propagation  of  a  diffuse 
interface  with  approximate  allowance  for  stress  variation  across  the  interface.  In  thermodynamic  equilibrium,  s2  =  0,  and 
Eq.  (36)  reduces  to  er  =  erM  and  t  =  xM. 

The  dependence  of  erM  -  a  on  aM  for  several  values  of  tm  is  presented  in  Fig.  3.  Despite  significant  differences  between  the 
2-3-4  and  2-4-6  potentials,  the  <rM  -  a  versus  <rM  curves  for  the  two  potentials  nearly  coincide.  The  corresponding  relation 
between  interface  velocity  and  aM  for  the  2-3-4  potential,  determined  with  the  help  of  Eq.  (13)  (for  a  =  <rM)  and  (37)  is 
shown  in  Fig.  4.  In  addition  to  the  parameter  values  for  NiAl  from  Eq.  (26),  we  also  used  the  values  y  =  0.2,  g  =  0.1,  and 
9  =  300  K  for  Figs.  3  and  4. 


4.  Numerical  study:  propagation  and  stability  of  A-M,  M-M,  and  M-A-M  interfaces 


In  this  section  we  present  the  results  of  FEM  simulations  of  A-M  and  M-M  interface  propagation,  interactions  between 
interfaces,  and  the  stability  of  interface  propagation.  Introducing  the  dimensionless  parameters 


Ccl  —  Gcl/A.  x  —  x\j~p'  t  — 

where  A  =  A0(9e  -  0C)  (defined  in  Eq.  (18)),  we  obtain 

dri  SGg.  (dC  d2?A 

dt  Sr]  \ydr]  dx2)' 


(39) 


(40) 


We  use  boundary  conditions  at  the  ends  of  the  sample  corresponding  to  equal  surface  energies  of  A  and  M: 

dgtlL 2)  dfim)  0  (41) 

dx  dx  ’  ' 

where  l  is  the  sample  length.  The  material  parameters  are  those  of  NiAl  given  in  Eq.  (26);  in  addition,  we  have  used  a  =  2.98. 
The  temperature  is  fixed  at  9  =  9e;  the  driving  force  for  phase  transformation  is  varied  by  changing  the  stress  a. 


Fig.  3.  Dependence  of  aM  -  a  on  oM  for  four  values  of  the  shear  stress  z  for  the  2-3-4  potential:  (1)  z  =  0;  (2)  z  =  100  MPa;  (3)  z  =  500  MPa;  (4) 
1  =  1000  MPa. 
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Fig.  4.  Relation  between  interface  velocity  and  aM  for  the  2-3-4  potential:  (1)  t  =  0;  (2)  t  =  100  MPa;  (3)  t  =  500  MPa;  (4)  r  =  1000  MPa. 


For  temperatures  around  the  phase  equilibrium  temperature  0e,  an  A  region  appears  between  martensitic  variants;  we 
study  this  case  as  well.  When  stresses  reach  and  exceed  the  value  corresponding  to  lattice  instability,  we  expect  that  homo¬ 
geneous  or  surface-induced  nucleation  occurs  in  addition  to  interface  propagation,  and  it  may  interact  with  the  interface 
propagation. 

4.3.  A- M  interfaces 

In  Fig.  5a  and  b  we  show  the  evolution  of  two  layers  of  M  separated  by  A  from  a  stationary  state  at  a  =  0,  under,  respec¬ 
tively,  tensile  and  compressive  stresses  of  magnitude  500  MPa.  In  tension  (si  =  1 .431  GPa,  s2  =  -1 .290  GPa,  thus  metastable 
A  and  stable  M)  the  M  layers  broaden  by  interface  propagation  until  the  A  has  completely  disappeared,  while  in  compression 
(si  =  2.072  GPa  and  s2  =  1.290  GPa,  metastable  M  and  stable  A)  the  M  layers  narrow,  eventually  disappearing.  Note  that  A  is 
unstable  (si  =  0)  for  er  =  2.733  GPa. 

In  the  next  simulation,  the  initial  condition  is  again  an  A-M  diffuse  interface  in  equilibrium  at  a  =  0.  At  a  =  2.733  GPa 
and  even  at  a  =  2.9335  GPa,  where  A  is  unstable,  the  interface  propagates  into  the  A  region  (si  <  0,  A  — >  M),  but  homoge¬ 
neous  transformation  of  the  unstable  A  does  not  occur.  This  is  because  8G/dr\  =  0  for  )/  =  0  in  our  model  so  there  is  no  local 
driving  force  for  the  transformation  of  A.  To  initiate  the  instability  we  impose  the  perturbation  ij(x)  =  0.01  sin(80nx).  For 
a  =  2.5335  GPa  (si  >  0,  s2  <  0.  A  metastable)  the  perturbation  rapidly  disappears,  and  for  a  =  2.9335  GPa  (Si  <  0,  A  unsta¬ 
ble)  the  initial  perturbation  again  disappears  except  near  the  A  surface  of  the  specimen  where  a  M  nucleus  appears  and 
grows  (the  nucleus  energy  at  the  surface  is  half  that  in  the  bulk  -  Levitas  et  al.,  2003);  see  Fig.  6a.  After  complete  transfor¬ 
mation  to  M  at  the  surface,  the  newly  formed  A-M  interface  propagates  toward  the  first  (Fig.  6a).  After  the  interfaces  meet  at 
the  point  ij  =  0,  they  continue  to  propagate  into  one  another  until  the  A  completely  disappears.  The  same  process  occurs  at 
the  limit  of  A  instability  (a  =  2.733  GPa),  but  the  newly  formed  interface  moves  at  a  lower  speed. 

If  the  initial  perturbation  is  of  the  form  ;;  =  0.01  |sin(407tx)|  (non-symmetric  with  respect  to  /;  =  0)  then  at  a  =  2.5335  GPa 
(A  metastable),  in  addition  to  the  propagating  interface,  the  A  homogeneously  transforms  to  M.  At  a  =  2.733  GPa  (si  =  0) 
and  a  =  2.9335  GPa  (si  <  0,  A  unstable)  the  homogeneous  transformation  occurs  faster  (Fig.  6b).  If  the  sample  is  sufficiently 
long  then  the  A  disappears  (homogeneously)  before  the  interface  reaches  the  surface. 


x/1  xll 


Fig.  5.  Evolution  of  austenite-martensite  microstructures.  The  initial  condition  is  a  stationary  state  of  two  layers  of  martensite  at  o  =  0  for  l  =  27.3  nm.  (a) 
Under  an  applied  stress  a  =  500  MPa  the  layers  broaden  by  motion  of  all  four  interfaces  until  the  transformation  is  complete,  (b)  For  o  =  -500  MPa,  the  M 
layers  narrow  and  finally  disappear. 
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Fig.  6.  Instabilities  in  the  propagation  of  austenite-martensite  interfaces.  The  initial  condition  is  a  stationary  austenite-martensite  microstructure  at  a  =  0 
for  /  =  27.3  nm.  (a)  The  perturbation  ii(x)  =  0.01  sin(80nx)  is  applied.  For  a  =  2.9335  GPa  (stress  exceeding  the  lattice  stability  limit),  in  addition  to  interface 
propagation,  a  martensitic  nucleus  appears  and  grows  at  the  sample  surface.  After  complete  transformation  to  M  at  the  surface,  which  forms  a  second  A-M 
interface,  the  two  interfaces  move  toward  one  another  until  the  martensitic  transformation  is  complete,  (b)  Application  of  the  perturbation 
Y][x)  =  0.01|sin(407rx)|.  For  a  =  2.9335  GPa  (A  unstable),  along  with  interface  propagation,  the  austenite  homogeneously  transforms  to  martensite. 


4.2.  M_— M+  and  M_— A— M+  interfaces 

In  Fig.  7  we  show  a  twinned  M_  — M+  microstructure,  stable  at  r  =  0.  After  applying  the  shear  stress  z  =  500  MPa  the 
microstructure  evolves  via  interfaces  propagation  into  a  homogeneous  M+  final  state.  Note  that  the  normal  stress  a  does 
not  contribute  to  the  driving  force  for  M_— M+  propagation.  The  single  interfaces  propagate  for  any  applied  stress. 

As  A-M  equilibrium  (AG"  =  0)  is  approached,  an  M_— M+  interface  splits  into  M_— A  and  A— M+  diffuse  interfaces  sepa¬ 
rated  by  a  layer  of  A  (soliton  splitting  -  Falk,  1983;  Levitas  et  al.,  2003).  In  Fig.  8  we  display  the  results  of  our  FEM  simulations 
of  a  propagating  M_— A— M+  interface  for  8  =  0e.  We  assume  y  =  0.2,  hence  M_  loses  its  stability  at  z  =  A/(6-  a)y  = 
2.899  GPa,  and  A  loses  its  stability  at  z  =  A/ay  =  2.938  GPa,  very  close  values.  Initially,  z  =  0  and  the  interface  is  stationary, 
but  under  the  applied  shear  stress  z  =  2.689  GPa  (M_,  A  metastable)  the  interface  propagates  and  the  width  of  the  A  layer 
concurrently  widens  (Fig.  8a).  The  reason  for  this  is  that  the  M_— A  interface  (metastable-metastable)  propagates  faster  than 
the  A— M+  interface  (metastable-stable);  the  propagating  M_— A— M+  interface  is  stable  (no  A  growth)  if  all  phases  are  stable 
or  metastable. 

To  check  for  an  instability  in  the  interface  propagation  due  to  A  instability,  an  initial  perturbation  rj(x)  =  0.04|sin  607tx| 
was  used  in  the  A  layer.  At  z  =  2.689  GPa  (M_,  A  metastable),  the  perturbation  disappears,  interface  propagation  is  stable, 
and  the  width  of  the  A  layer  increases,  as  in  Fig.  8a.  Under  z  =  2.938  GPa  (M_  unstable,  lower  limit  of  A  instability),  the  sit¬ 
uation  is  very  similar,  but  there  are  small  deviations  of  r;  from  zero  in  the  A  layer  (Fig.  8b).  For  z  =  3.1885  GPa  (A,  M_  unsta¬ 
ble)  a  section  of  the  A  layer  transforms  to  M+  during  interface  propagation.  This  transformation  runs  to  completion,  thus 
forming  a  new  stable  A— M+  interface  that  continues  to  propagate;  shown  in  Fig.  8c. 

The  stability  of  the  variant  M_  under  varying  stress  was  studied  for  the  initial  perturbation  r](x)  =  0.04|sin607tx|  in  the  M_ 
region.  At  z  =  2.689  GPa  (M_,  A  metastable)  the  perturbation  had  no  effect.  At  t  =  2.938  GPa  (loss  of  A  stability)  and 
t  =  3.1885  GPa  (Fig.  8d),  the  M_  homogeneously  transforms  to  A  as  the  A— M+  interface  propagates. 


x!l 


Fig.  7.  The  initial  state  is  a  stationary  twinned  martensitic  microstructure  at  z  =  0,  /  =  13.6  nm.  For  an  applied  stress  of  t  =  500  MPa,  the  twins  disappear, 
resulting  in  a  homogeneous  M+  final  state. 
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Fig.  8.  Propagation  of  a  M_— A— M+  interface  and  its  instability;  /  =  20.4  nm.  (a)  For  the  shear  stress  t  =  2.689  GPa  (M_,  A  metastable),  the  width  of  the  A 
region  widens  as  the  interface  propagates,  (b)  For  the  shear  stress  t  =  2.938  GPa  (M_  unstable,  minimum  stress  for  A  instability),  the  initial  perturbation 
rj(x)  =  0.04|sin(607tx)|  in  the  A  region  disappears  and  the  interface  propagation  is  stable,  (c)  For  the  shear  stress  x  =  3.1885  GPa  and  the  initial  perturbation 
rj(x)  =  0.04|sin607rx|  in  A,  a  portion  of  the  A  transforms  to  M+  during  interface  propagation.  This  transformation  results  in  a  new  stable  A— M+  propagating 
interface,  (d)  For  the  shear  stress  t  =  3.1885  GPa  and  the  initial  perturbation  r\(x)  =  0.04|s/n(607cx)|  in  the  M_  region,  the  M_  homogeneously  transforms  to 
A  as  the  A— M+  interface  propagates. 


5.  Three-dimensional  Landau  model  for  multivariant  stress-induced  martensitic  phase  transformations  and  2D 
simulations 

Recently  (Levitas  and  Preston,  2002a, b;  Levitas  et  al.,  2003),  a  3D  Landau  theory  for  stress-induced  martensitic  phase 
transformations  was  developed.  It  can  describe  phase  transformations  between  A  and  M,  and  between  M,  variants  for  any 
crystal  structures,  as  well  as  typical  stress-strain  curves  for  phase  transformations  in  shape  memory  alloy,  steels,  and  inter- 
metallics.  Also,  it  can  incorporate  all  temperature-dependent  thermomechanical  properties  of  the  A  and  M  phases.  The  Gibbs 
energy  G  and  the  transformation  strain  tensor  £t  are  of  the  form  (Levitas  and  Preston,  2002a, b) 

1  n  n- 1  n 

C=--<r:l:<r-a:£t-a:£d  +  ^/(0,  rjk)  +  Y.  *lj), 

k=  1  i=l  j=i+l 

=  +  (4  -  2a^l  +  (a  -  -  Y  Y  +  rijLj,), 

k=  1  i=l  j=i+l 

with 

nit)  =  H  +  (4AC"  -  2A)ri  +  (A  -  3AG0)rjl  Uj  =  (a  -  3)4  +  38* , 

Fijilr  >1j)  =  WjO  ~  >h  ~  -  ilj)2  ~  >h  ~  >ljl  +  Dih'lj}  + 

Here,  a  is  the  stress  tensor;  17,  and  4  are  the  order  parameter  and  the  transformation  strain  tensor  of  the  ith  variant,  where 
i  =  0  corresponds  to  A  and  =  0;  k  is  the  elastic  compliance  tensor  (assumed,  for  simplicity,  to  be  the  same  for  A  and  M,); 
and  A,  A.  B,  and  C  are  material  parameters.  In  comparison  with  equations  in  Levitas  and  Preston  (2002a, b)  and  Levitas  et  al. 
(2003),  we  added  the  eigen  strain  £d  due  to  crystal  defects  (dislocation,  points  defects,  grain  and  subgrain  boundaries)  that 
we  will  need  in  the  next  section.  The  parameters  A  and  A  characterize  the  thresholds  for  A  +-►  M,  and  M,-  <->  M,  transforma¬ 
tions,  while  B  and  C  control  the  Gibbs  energy  away  from  both  the  A  and  M,  minima  and  the  minimum-energy  paths  between 


and 

+  rij)(A-A). 


(42) 

(43) 


(44) 
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the  minima,  therefore  they  do  not  affect  the  phase  equilibrium  and  transformation  conditions.  The  constants  B  and  D  do  not 
contribute  to  phase  equilibrium  and  instability  conditions.  They  are  used  to  avoid  non-physical  energy  minima  that  may  ap¬ 
pear  for  such  a  complex  polynomial. 

The  evolution  of  the  n  order  parameters  is  described  by  the  n  coupled  TDGL  equations 


& U 

dt 


(45) 


where  Ggl  =  fyCadQ  and  Ga  =  G  +  ( Vi/fc)2  is  the  GL  energy.  In  the  following  we  use  Eq.  (16)  for  AG"  and  A,  the  mate¬ 

rial  parameters  for  the  cubic-tetragonal  phase  transformation  in  NiAl  found  in  Levitas  and  Preston  (2002a, b)  and  Levitas  et  al. 
(2003)  and  summarized  in  Eqs.  (26)  and  (28),  as  well  as 


a  =  2.980,  A  =  5320  MPa,  B  =  0,  D  =  5500  MPa,  v  =  0.238,  (46) 

where  v  is  Poisson’s  ratio.  For  our  2D  FEM  simulations  we  consider  just  two  of  the  three  possible  NiAl  martensitic  variants 
e]  =  (0.215;  -0.078;  -0.078),  s2t  =  (-0.078;  0.215; -0.078),  (47) 

The  potential  (42)  leads  to  the  elasticity  relation 

E  =  7,  :  G  +  Et  -h  Erf.  (48) 

We  solve  Eq.  (48)  for  the  stress  tensor 

a  =  k  1  :  (£  —  Et)  +  ad,  (49) 

where  we  designated  ad  =  -2T1  :  £d  for  an  oscillatory  defect  stress  field,  which  we  will  discuss  in  detail  in  the  next  Section.  In 
addition,  the  standard  equilibrium  equations  and  the  relationship  between  strains  and  displacements  u  are  given  by: 

V  •  a  =  0,  £  =  (Vu)j,  (50) 

where  (••■)s  means  symmetrization.  Thus,  the  complete  system  of  equations  describing  the  phase  transformation  is  com¬ 
prised  of  Eqs.  (42)-(45)  and  (49),  (50).  Eqs.  (43),  (45),  (49),  and  (50)  are  similar  in  structure  to  the  coupled  equations  of  ther¬ 
moelasticity.  Eq.  (45)  resembles  a  set  of  heat  conduction  equations  for  n  temperatures  with  temperature-dependent  heat 
sources  XdG/dr\k.  The  transformation  strain  £t  corresponds  to  a  thermal  strain  with  a  complicated  temperature  dependence 
(43).  This  correspondence  between  phase  transformation  and  thermoelasticity  equations  has  important  computational  con¬ 
sequences:  finite  element  thermoelasticity  codes  can  be  used,  after  some  minor  modifications,  for  phase-field  model  simu¬ 
lations  of  phase  transformations.  In  contrast  to  approaches  based  on  the  spectral  (fast  Fourier  transform)  method  (Artemev 
et  al.,  2001;  Curnoe  and  Jacobs,  2001  a, b;  Jacobs  et  al.,  2003;  Jin  et  al.,  2001 ;  Lookman  et  al.,  2003a, b;  Rasmussen  et  al.,  2001; 
Seol  et  al.,  2003;  Wang  and  Khachaturyan,  1997;  Wang  et  al.,  2001),  the  finite  element  approach  allows  us  to  easily  expand 
the  treatment  to  heterogeneous  materials,  large  strains,  arbitrary  boundary  conditions,  and  complex  material  models. 
Because  the  potential  (42),  (43)  accurately  describes  the  important  features  of  martensitic  phase  transformations,  we  expect 
that  our  calculated  microstructure  evolution  in  this  study  is  more  realistic  than  that  predicted  by  other  approaches  (Artemev 
et  al.,  2001;  Curnoe  and  Jacobs,  2001  a, b;  Jacobs  et  al.,  2003;  Jin  et  al.,  2001 ;  Lookman  et  al.,  2003a,b;  Rasmussen  et  al.,  2001; 
Seol  et  al.,  2003;Shenoy  et  al.,  1999;  Wang  and  Khachaturyan,  1997;  Wang  et  al.,  2001).  In  this  work  we  used  the  code  FEAP 
(Zienkiewicz  and  Taylor,  2000). 

We  consider  a  square  sample  of  size  /  =  60  nm  under  plane  strain  conditions.  The  temperature  is  0  =  100  K,  and  at  the 
boundaries  constant  stresses  ab  =  (eq,, ,  ob2)  are  applied  which  satisfy  n  ■  Vt/k  =  0,  where  n  is  the  normal  to  the  boundary. 
In  our  simulations,  10,000  eight-node  quadrilateral  finite  elements  were  used.  The  thermodynamic  driving  force  for  vari¬ 
ant-variant  phase  transformations  due  to  the  boundary  conditions  is  X"^2  =  ab :  (e2  -  ej)  =  (ab2  -  (Tbr)(£22  -  £^)  = 
0.293(<Ti,2  -  or i).  The  thermodynamic  driving  force  for  the  A  — >  M,  phase  transformation  due  to  boundary  conditions  is 
Xj_M-  =  ob  :  £*  =  OttE'n  +  ob2£\2.  Note  that  large  internal  stresses  that  exceed  the  stability  limit  of  austenite  may  lead  to  val¬ 
ues  17,  >  1  where  the  Landau  potential  may  exhibit  unphysical  behavior  (Levitas  and  Preston,  2002a, b;  Levitas  et  al.,  2003).  To 
avoid  this,  we  enforced  the  constraint ^  1  in  our  numerical  simulations. 

In  our  first  simulation,  the  initial  conditions  are  randomly  distributed  and  r;2  fields  (Fig.  9),  and  the  boundary  condi¬ 
tions  are  abi  =  10  GPa  (an  arbitrary  stress  that  provides  a  large  driving  force  for  the  phase  transformation  A  — >  M,);  therefore 
Xj^2  =  0.  The  solutions  of  Eqs.  (42)-(45)  and  (49),  (50)  are  shown  in  Fig.  9:  the  first  column  shows  the  evolution  of  ijv  the 
second  and  third  columns  show  the  stresses  a ,  and  oq  -  er2  (since  the  thermodynamic  driving  force  X2^i  is  proportional  to 
o" i  -  a2),  and  the  last  column  depicts  the  local  driving  force,  dG/dij ,,  for  the  evolution  of  jj,  .  The  evolution  of  ij2  and  dG/dij2 
are  qualitatively  the  same;  for  example,  in  the  bottom  figure  for  >/, ,  the  dark  blue  field  corresponds  to  r/2  =  1,  i.e.,  with  the 
variant  M2.  After  passing  through  a  complex  microstructure  containing  plates,  laths,  and  needles  the  solution  converges  to  a 
primitive  M,— M2  twin  microstructure  (Fig.  9).  Stresses  and  local  driving  forces  are  concentrated  at  the  interfaces,  and  local 
driving  forces  are  quite  low.  However,  this  microstructure  is  not  stable  -  driven  by  a  decrease  in  gradient  energy  through  a 
reduction  in  interface  area  it  converges  to  a  single  phase,  M,  or  M2,  depending  upon  initial  conditions.  The  propagation  direc¬ 
tion  of  each  planar  interface  as  well  as  the  final  state  variant  can  be  easily  predicted:  each  interface  moves  to  reduce  its 
length,  thus  the  variant  occupying  the  central  part  of  the  sample  is  the  survivor. 
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Fig.  9.  Evolution  of  initially  randomly  distributed  fields  and  r/2  (first  column),  rj,  and  -  a2  (second  and  third  columns),  and  the  local  driving  force, 
aC/Sr/,,  for  evolution  of  /j1  (right  column)  for  rrM  =  ab2  -  10  GPa.  The  microstructure  finally  converges  to  homogeneous  M2  with  zero  internal  stresses  and 
driving  forces. 


Note  that  the  normals  to  the  interfaces  vary  from  parallel  to  the  boundary  in  a  very  thin  surface  layer  (due  to  the  bound¬ 
ary  condition  n  ■  Vi;fc  =  0)  to  45°  to  the  boundaries  (determined  by  elastic  energy  minimization).  This  is  typical  of  all  sim¬ 
ulations  below  but  difficult  to  see  in  the  figures. 

In  our  second  simulation,  the  initial  conditions  are  =  iy2  =  0.1  in  a  circle  of  radius  2  nm  (an  embryo)  at  the  center  of  the 
sample  and  zero  elsewhere.  The  boundary  stresses  are  erM  =  ab2  =  15  GPa;  hence,  the  system  is  initially  symmetric  with  re¬ 
spect  to  each  variant.  A  very  complicated  microstructure  consisting  of  a  combination  of  four  herring  bone  types  of  micro¬ 
structure  is  developed  (Fig.  10).  At  some  stage,  the  equivalence  of  both  variants  is  violated-units  of  M2  coalesce  at  the 
center,  and  Mi  coalesces  in  four  triangles  near  the  sample  boundaries.  Other  units  evolve  to  a  needle-like  microstructure 
(similar  to  that  observed  experimentally  in  Abeyaratne  et  al.  (1996)).  Later,  the  solution  evolves  into  a  primitive  twinned 
microstructure  and  finally  to  homogeneous  M2. 

Thus,  as  discussed  in  the  Introduction,  in  both  of  these  simulations,  for  which  the  boundary  stresses  are  specified,  the  lack 
of  athermal  friction  leads  asymptotically  to  a  single  phase. 


6.  Athermal  threshold  for  interface  propagation 

The  above  equations  are  valid  when  the  magnitude  of  the  driving  force  exceeds  all  barriers  due  to  point  and  line  defects, 
as  well  as  the  Peierls  barrier,  because  these  barriers  were  not  included  in  the  preceding  analysis.  At  these  high  driving  forces 
the  interface  motion  is  governed  by  the  phonon  (and  at  lower  temperatures,  by  the  electron)  drag  mechanism,  but  the  bar¬ 
riers,  which  result  in  an  athermal  friction  force  opposing  the  interface  motion  (Ghosh  and  Olson,  1994;  Grujicic  et  al.,  1985; 
Olson  and  Cohen,  1986),  must  be  taken  into  account  at  lower  driving  forces.  As  discussed  in  the  Introduction,  the  inclusion  of 
an  athermal  threshold  in  the  phase  field  approach  is  not  straightforward  because  the  phase  field  approach  does  not  involve  a 
separate  equation  for  the  interface  but  rather  an  evolution  equation  for  the  order  parameter.  We  can  put,  for  example,  2  =  0 
if  |S2/12|  <  I<,  however,  this  can  arrest  unphysical  intermediate  configurations,  thus  the  system  never  converges  to  a  realistic 
microstructure  consisting  of  A  and  M  variants  separated  by  diffuse  interfaces,  either  stationary  or  moving.  We  now  propose 
two  schemes  for  solving  this  problem. 
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Fig.  10.  Evolution  of  an  embryo,  rjt  =  0.1,  in  a  circle  of  radius  of  2  nm  at  the  center  of  the  sample  for  er*,,  =  15  GPa.  Left  column:  y\x  ;  second  and  third 
columns:  ox  and  Ci  -  o2  \  right  column:  dG/drj^,  the  local  driving  force  for  evolution  of  rj^.  The  stationary  microstructure  is  M2.  In  the  first  two  rows,  red  is 
Mi,  and  the  blue  outside  of  the  Mi  regions  is  A;  the  distribution  for  M2  is  rotated  90°  from  the  distribution  for  Mi.  In  the  last  four  rows,  red  is  Mi  and  blue  is 
M2.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


6.1.  Athermal  threshold  I:  modified  evolution  equations 


To  introduce  a  threshold  for  interface  propagation,  we  need  to  find  a  way  to  distinguish  between  interface  propagation 
and  an  arbitrary  non-stationary  process.  At  all  points  of  the  propagating  interface  (drj  /  dt)  /  (drj  /  dx)  =  ±c;  therefore,  the  spa¬ 
tial  derivative  of  this  expression  is  zero,  which  results  in 


d2r\  dt]  drj  cfrj  _ 
dxdt  dx  dt  dx2  ~ 


(51) 


Near  the  interface,  |K|  <c  1 .  We  can  account  for  interface  friction  by  requiring  very  slow  dynamics,  which  is  realized  for 
very  small  values  of  the  kinetic  coefficient  X,  in  narrow  bands  around  the  interfaces  if  the  driving  force  is  less  than  a  critical 
value,  /(: 
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X  =  o  for  [J?|  <  g  and  [s2/12|  <  K,  X  =  X0  otherwise. 


(52) 


Here  q  is  a  small  quantity  of  order  the  computational  error,  and  o  is  sufficiently  small  that  dr\/dt  is  negligible.  We  use  o  in¬ 
stead  of  zero  to  avoid  X  =  0  at  points  where  dr]/dt  =  drj/dx  =  0  (R  =  0),  but  which  do  not  belong  to  the  moving  interface,  e.g., 
at  the  maximum  of  a  critical  M  nucleus. 

Fig.  11  shows  the  results  of  a  numerical  simulation  of  the  evolution  of  an  initial  sinusoidal  order  parameter  (2-4-6  po¬ 
tential,  AC0  =  0).  First,  two  M_(;/  =  — 1)— M+(i)  =  1)  interfaces  are  formed  but  their  motion  is  arrested  because  |s2/l 2|  <  K. 
The  magnitude  of  the  driving  force  is  then  increased  above  the  threshold  K  (by  increasing  r)  and  the  interfaces  begin  moving 
toward  each  other.  They  first  meet  at  the  point  (x,  t /)  =  (0,  -1).  Subsequently, )?( 0,  t)  increases  toward  )/  =  +1  until  M_  com¬ 
pletely  disappears,  resulting  in  a  final  state  of  homogeneous  M+.  It  is  clear  that  our  goal  has  been  accomplished:  interface 
motion  can  be  both  arrested  and  released,  while  intermediate  unstable  configurations  are  not  affected  by  the  athermal 
threshold. 

Eq.  (52)  has  the  drawback  that  the  kinetic  coefficient  1  and  the  interface  velocity  jump  from  zero  to  a  finite  value  at 
|s2/l 2|  =  K.  However,  experimental  observations  (Ghosh  and  Olson,  1994;  Grujicic  et  al„  1985;  Levitas  et  al.,  2002c;  Olson 
and  Cohen,  1986)  show,  and  sharp  interface  theory  predicts,  that  the  magnitude  of  the  interface  velocity  depends  on  the 
excess  of  the  driving  force  over  the  athermal  threshold,  s2/12  -  I<  for  s2/12  >  I<  and  s2/12  +  K  for  s2/12  <  —K,  rather  than 
on  [s2/12|.  Note  that  these  replacements  are  equivalent  to  a  ■.  £t  — >  a  :  £t  ±K.  If  the  athermal  threshold  is  caused  by  the  stress 
fields  of  various  crystal  defects,  then  all  occurrences  of  a  :  et  in  Eqs.  (13)  and  (14)  for  c4  and  c6  must  be  replaced  by  a  :  e,  ±  K. 
Thus,  the  interface  speeds  for  |s2/12|  >  I<  are 


c4  =  Xn^Jp 


s2(1-JC/|s2/12|) 


Cb  =  Xg  \ 


\/4s i  -  s2  +  sign(s2) 4(3  -  a)K 
Tp  s2(\~K/\s2/U\) 

3  ^3s,  -  s2  +  sign{s2) 3(4  -  a)K 


(53) 

(54) 


The  interface  speeds  Eqs.  (53)  and  (54)  can  be  incorporated  in  our  phase  field  model  through  X0  (kinetic  coefficient  for 
|s2/l 2|  >  K\  see  Eq.  (52))  defined  by  requiring  that  c4  and  c6  as  given  by  Eqs.  (13)  and  (14)  with  1  X0  are  equal  to  the  cor¬ 
responding  velocities  in  Eqs.  (53)  and  (54): 


A  =i  1  ~/C/|s2/12| 

04  4  v/1  +  4sign(s2)(3  -  a)K/{Asx  -  s2) 

X  =;.  1  -7C/|s2/12] 

06  6  v/1  +  3sign(s2)(4  -  a)/C/(3s!  -  s2) 


(55) 

(56) 


For  the  general  case  of  three  dimensions  and  multivariant  phase  transformations  described  by  m  order  parameters  f/,, 
(9r/i/9t)/|V)/l|  =  ±c  on  the  interface.  The  local  direction  of  interface  propagation  is  n  :=  S7r/i/\S7r]j\‘,  let  r  be  the  local  coordi¬ 
nate  along  n.  Then  (drij dt) / {di^dr)  =  ±c  at  all  points  of  the  interface.  Again  setting  the  r  derivative  to  zero  yields 


d\dr]1_dr]1cftj1  = 
drdt  dr  dt  dr 2 


(57) 


For  each  i  we  can  express  Xt  in  the  same  way  as  in  Eq.  (52). 


x// 

Fig.  11.  Numerical  simulation  of  the  evolution  of  the  initial  configuration  (curve  1 )  rj(x)  =  cos(7tx/2)  for  the  threshold-type  kinetics  defined  by  Eq.  (52).  The 
initial  configuration  evolves  into  two  M_— M+  interfaces  (curve  2)  which  are  arrested  due  to  the  condition  |s2/12|  <  K.  Subsequently,  |s2/12|  is  increased 
above  the  threshold  K  and  interfaces  move  toward  each  other.  They  first  meet  at  the  point  (x,  rj )  =  (0,  -1)  and  then  77(0,  t)  increases  toward  rj  =  +1  leading 
to  complete  disappearance  of  M_  and  to  a  homogeneous  M+  structure. 
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Our  approach  can  also  be  used  for  thermally  activated  interface  motion,  which  is  described,  for  example,  by  the  equation 
(Ghosh  and  Olson,  1994;  Grujicic  et  al.,  1985) 


c  =  cs  exp 


^(\-[(\s2/U\-K)/W0f)b 


(58) 


where  Q0  is  the  activation  energy,  W0  is  the  height  of  the  driving  force  barrier  corresponding  to  the  rate-controlling 
(short-range)  obstacles,  k  is  the  Boltzmann  constant,  and  p  and  b  are  constants.  All  of  the  material  parameters  appearing 
in  (58)  have  been  estimated  for  Fe-22.31  Ni-2.888Mn  alloy  (Table  1  in  Levitas  et  al.,  2002c).  As  above,  we  can  define 
threshold-dependent  kinetic  coefficients  X04  and  X06  by  equating  c4  and  c6  as  given  by  Eqs.  (13)  and  (14)  with  X^  X0  to 
c  from  Eq.  (58). 

Note  that  for  AG"  -  ere  =  0,  the  limit  condition  for  the  A  — >  M+  transformation,  XA^M+  =  — s2/12  =  zy  =  I<,  and  that  for  the 
M_  — >  A  transformation,  XM _^a  =  s2/12  =  zy  =  K,  result  in  Xm_^m+ =  XA^M+ +XM_^A  =  2zy  =  2I<,  i.e.,  the  threshold  for 
M_  — >  M+  interface  propagation  is  2 K.  In  this  case  a  finite  A  region  appears  between  the  martensitic  variants,  and  there  is 
no  need  to  introduce  a  separate  threshold  condition  for  variant-variant  transformations.  We  consider  only  the  case  where 
the  interface  moves  into  M_;  the  opposite  case  leads  to  identical  conditions  for  interface  propagation.  If  AG"  -  cre^O,  then 
XA  m^XXm  ,A  and  it  is  possible  that,  for  example,  XM_^A  <  0,  in  which  case  the  transformation  condition  XM_^A  >  I<  for 
the  M_  -  A  interface  cannot  be  satisfied.  This  means  that  one  must  consider  the  M_  — >  M+  interface  separately  and  apply 
the  condition 

2  =  o  for  |i?|  <  q  and  |ry|  <  I(,  X  =  X0  otherwise.  (59) 


6.2.  Athermal  threshold  II:  oscillatory  defect  stress  fields 

Usually,  kinetic  equations  for  thermally  activated  dislocation  (Kocks  et  al.,  1975)  or  interface  motion  are  derived  by 
considering  spatially  oscillatory  energy  profiles  that  represent  the  Peierls  barrier  or  the  energy  profiles  of  various  defects 
(Ghosh  and  Olson,  1994).  We  will  introduce  a  spatially  oscillatory  energy  profile  through  oscillatory  stress  fields,  ad,  due 
to  defects;  the  stress  tensor  is  given  by  Eq.  (49).  If,  for  example,  ad  =  Ken  sin(2nx/b)/(Bn  ■  £ti),  then  the  Landau  potential  con¬ 
tains  an  additional  transformation  work,  Ksin{2nx/b),  which  produces  an  effect  similar  to  an  interface  friction  l<  for  motion 
in  both  directions.  A  similar  effect  can  be  obtained  by  adding  Ksin(2nx/b)  to  the  change  in  the  thermal  part  of  the  Gibbs 
energy,  AG". 

The  period  b  is  determined  by  the  actual  field  of  defects.  If  the  origin  of  the  athermal  threshold  is  the  Peierls  barrier,  then  b 
is  the  period  of  the  crystal  lattice.  If  the  stress  field  is  due  to  stochastically  distributed  dislocations,  then  the  period  is  the 
averaged  distance  between  dislocations,  a  well-known  length  scale.  Multiple  length  scales  would  be  present  in  a  more  de¬ 
tailed  treatment  of  the  system  of  dislocations. 

The  evolution  of  an  initial  sinusoidal  tj(x)  (2-4-6  potential)  under  the  shear  stress  z  =  75  MPa  is  shown  in  Fig.  12.  In  the 
absence  of  the  periodic  field,  the  solution  evolves  to  M+,  but  in  contrast,  a  M_— A— M+  microstructure  (Fig.  12a,  configuration 
3)  is  formed  in  the  presence  of  the  periodic  stress  field  rd  =  ~I00sin(8nx)  MPa  (Fig.  12c).  After  increasing  the  stress  z  to 
125  MPa,  the  M_— A—  M+  interface  propagates  to  the  left  (Fig.  12b). 

7.  Microstructure  evolution  in  the  presence  of  oscillatory  defect  stress  fields 

In  all  calculations  we  use  the  material  parameters  for  the  cubic-tetragonal  transformation  in  NiAl.  We  perform  2D  sim¬ 
ulations,  impose  the  boundary  stress  Om  =  crb2  =  10  GPa,  take  0  =  100  K,  and  consider  only  the  two  NiAl  martensitic  variants 
with  transformation  strains  eJ  =  (0.215; -0.078; -0.078)  and  E2t  =  (-0.078;  0.215; -0.078).  For  ad  =  (<7dl,crd2),  the  thermo¬ 
dynamic  driving  force  for  the  Mi  — >  M2  phase  transformation  due  to  the  defect  stress  fields  is  X^_,2  =  ad  :  (s2  -  b})  = 
AEt((7d2  -  Gin),  where  Aet  =  e^2  -  e)2  =  e),  -  e2n  =  0.293.  The  thermodynamic  driving  force  for  the  A  — >  M,  phase  transforma¬ 
tions  due  to  the  stress  fields  of  the  defects  is  XA  >M  =  od:E\  =  ad^e'n  +  (7d2£j2.  In  the  particular  case  od2  =  crdl  one  has 
X?  >2  =  0  and  X"  M.  =  erdl  (ej,  +  e'2). 

We  study  four  classes  of  problems.  Some  preliminary  results  can  be  found  in  Levitas  and  Lee  (2007). 

1.  Solution  to  the  same  problem  as  in  Fig.  9  (random  initial  distribution  of  A  and  M)  but  with  the  periodic  defect  stress 
field  ad2  =  -cdl  =  10sin(167rx)sin(167ty)  and  consequently  X^2  =  20Aetsin(167rx)sin(167iy);  see  Fig.  13.  In  the  absence  of 
surface  energy  and  internal  stresses  due  to  transformation  strain,  a  minimum  of  the  energy  would  correspond  to  multi-con¬ 
nected  microstructures  with  M!  (M2)  in  the  regions  withX"^2  <  0  (X^_,2  >  0).  This  is  essentially  the  case  in  the  first  row  of 
Fig.  13.  The  elastic  energy  due  to  evolving  transformation  strain  fields  (within  evolving  martensitic  regions)  promotes  the 
appearance  of  straight  plates  aligned  at  45°,  and  surface  energy  reduction  drives  coalescence.  In  this  stage,  elastic  and  surface 
energy  exceed  the  potential  barriers  due  to  the  X"  >2  field,  and  the  units  coalesce  in  two  orthogonal  groups  of  twinned  plates. 
One  group  eventually  predominates,  resulting  in  a  stationary  microstructure  comprised  of  equidistant  twinned  plates.  Our 
periodic  fields  do  not  arrest  any  intermediate  microstructures  (0  <  r]t  <  1 ),  nor  do  they  change  the  orientation  of  the  mar¬ 
tensitic  plates  which  is  determined  by  minimization  of  elastic  energy,  but  they  do  slightly  distort  the  otherwise  planar 
interfaces. 
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Fig.  12.  Evolution  of  an  initial  sinusoidal  i/(x)  (curve  1)  in  the  periodic  defect  stress  field  rd(x)  =  -100sin(87rx)  MPa  shown  in  (c).  (a)  Under  a  shear  stress 
t  =  75  MPa,  a  M_—  A— M+  split  interface  is  formed  (curve  3).  (b)  After  increasing  the  shear  stress  to  125  MPa,  the  M_—  A—  M+  interface  propagates  to  the  left. 


Fig.  14  shows  the  solution  to  the  same  problem  as  in  Fig.  8  but  with  equal  periodic  stress  fields  of  defects 
<7d2  =  Crfi  =  10sin(167ix)sin(167ty)  (and  consequently  X^s2  =  0).  Again,  elastic  energy  due  to  evolving  transformation  strain 
fields  (within  evolving  martensitic  regions)  promotes  the  appearance  of  straight  plates  aligned  at  45°,  and  the  surface  energy 
promotes  coalescence. 

The  field  XdM ,_,A  affects  the  initial  stages  of  formation  of  the  microstructure,  reducing  the  initial  M,  units  to  sizes  compa¬ 
rable  to  its  period.  In  this  stage,  the  elastic  and  surface  energies  exceed  the  potential  barriers  due  to  the  XjJ,._,A  field,  and  the 
units  coalesce.  Once  the  system  has  converged  to  a  twin  microstructure  topologically  similar  to  that  in  Fig.  9,  it  is  arrested 
because  the  elastic  and  surface  energies  cannot  overcome  the  barriers  due  to  the  XjJ,._,A  field.  The  widths  of  different  mar¬ 
tensitic  plates  vary  by  several  times,  and  for  most  plates  the  width  exceeds  the  half  period  of  the  field  XjJ,._A.  Consequently, 
small  islands  of  incompletely  transformed  A,  stabilized  by  the  XjJ,._A,  appear  within  broad  martensitic  plates,  as  can  be  seen 
in  Fig.  14.  If  the  boundary  stresses  abi,  and  consequently  XA^M  ,  are  decreased,  then  these  islands  will  serve  as  heterogeneous 
nucleation  sites  for  A.  Although  the  volume  fractions  of  Mi  and  M2  vary  with  the  initial  conditions,  the  volume  fraction  of  Mi 
is  invariably  larger  than  that  of  M2.  The  presence  of  both  Mi  and  M2  is  due  entirely  to  the  defect  stress  field-if  removed,  the 
final  microstructure  consists  of  Mi  or  M2  only,  depending  on  initial  conditions.  Even  though  Xj_2  =  0,  the  stress  field  of  the 
defects  nevertheless  forms  a  mixed  Mi— M2  microstructure  because  it  arrests  the  formation  and  propagation  of  A— M, 
interfaces. 

2.  For  the  fields  ad i  =  5cos(167tx)cos(167ty)  and  ad2  =  5sin(167tx)sin(167ty),  X?_2  =  -5A£tcos(167i(x  +  y)),  which  is  peri¬ 
odic  along  the  diagonal.  The  initial  microstructure  evolves  into  a  perfect,  alternating,  equally  sized  Mi— M2  microstructure 
oriented  at  45°  (Fig.  15).  The  widths  of  the  M,  plates  are  determined  by  the  period  of  the  oscillating  stress  field-there  is 
no  lateral  coalescence  of  M,  plates. 

Next  -  see  Fig.  16  -  the  stationary  microstructure  obtained  in  the  previous  simulation  (Fig.  14)  was  taken  as  an  initial 
condition,  with  crw  =  10  GPa.  The  defect  stress  field  is  removed  (X^_2  =  0),  thus  there  is  no  athermal  threshold  or  hysteresis. 
After  removal  of  ad ,•  the  Mi  plates  widen  at  the  expense  of  M2  until  the  transformation  to  Mi,  driven  by  the  reduction  of  sur¬ 
face  energy,  is  complete. 

When,  prior  to  phase  transformation  completion,  the  stress  ab i  at  the  boundary  was  decreased  to  9.9  GPa,  the  reverse 
phase  transformation  Mi  — >  M2  occurred  by  equal  widening  of  all  M2  plates  until  the  transformation  was  complete. 

In  Fig.  17  we  plot  interface  velocity  v  vs.  the  macroscopic  driving  force  Xj_,2  for  four  interfaces  from  the  microstructure  in 
Fig.  16.  The  small,  non-zero  velocity  forX*  2  =  0  is  caused  by  the  driving  force  due  to  the  reduction  of  surface  energy;  with 
increasing  X^_,2,  the  contribution  of  the  surface  energy  is  negligible.  All  of  these  points  can  be  approximated  by 
v  =  -0.005  +  0.0034Xi^2  with  very  small  v  forX*  >2  =  0.  For  an  interface  of  constant  length,  v  =  0  for  X,  ^2  =  0,  which  is  also 
confirmed  by  our  ID  analytical  solution.  Thus,  the  athermal  threshold  and  hysteresis  are  absent. 
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Fig.  13.  Solution  to  the  same  problem  as  in  Fig.  9  but  with  the  periodic  defect  stress  field  ad2  =  nd]  =  10sin(16ra)sin(16jiy).  Left  column: ;  second  and 
third  columns:  rj,  and  <7,  cr2 ;  right  column:  the  local  driving  force,  dG/3/y, ,  for  the  evolution  of  ij,  .  The  bottom  row  is  the  stationary  solution  (red  is  M,  and 
blue  is  M2).  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


For  the  same  initial  and  boundary  conditions,  the  defect  stress  field  erdl  =  ad2  =  10sin(87t(x  +y))  leads  to  Xd_2  =  0  and  to 
XdA  iM  a  periodic  function  along  the  diagonal.  The  initial  microstructures  coalesce  into  alternating  Mr— M2  plates  of  different 
widths  oriented  at  45°;  Fig.  18. 

3.  In  Fig.  19  we  simulate  the  evolution  of  a  stochastic  initial  distribution  of  Mr  and  M2  in  a  defect  stress  field  for  which 
<jd2  -  <7dt  is  a  2D  plane  square  wave  (in  general,  the  periodic  extension  of  the  function  C[0(r)  -  fl(— r)]  on  -b  <  r  sg  b,  where  C 
is  a  constant,  r  is  the  coordinate  along  the  normal  to  the  wave,  and  2b  is  the  wavelength).  We  superimpose  equal  period,  in- 
phase  square  waves  for  ad2  and  -<rdl,  each  with  an  amplitude  of  2  GPa,  a  wavelength  of  0.2/,  and  normals  inclined  45° 
(Fig.  19).  The  amplitude  of  the  ad2  -  <rdi  square  wave  is  4  GPa  (|Xf  ,2  =  Aet|<Td2  — ■  <rdl  |  =  4Asf  GPa  =  1.172  GPa),  thus  the  ex¬ 
pected  stress  hysteresis  is  H  =  \ od2  —  crdl |  =  4  GPa.  Indeed,  keeping  ob2  =  10  GPa  for  6  <  crM  <  14  GPa  the  interfaces  do  not 
move;  for  obi  =  10  GPa  the  system  converges  to  a  twinned  microstructure  with  twins  occupying  the  regions  of  positive  and 
negative  X^2.  Only  for  abl  >  14  GPa  (ffbl  <  6  GPa)  does  the  variant  Mi  (M2)  grow  until  the  phase  transformation  is  com¬ 
plete.  Thus  plane  square  wave  defect  stress  fields  produce  a  finite,  rate-independent  athermal  hysteresis  (deviation  from 
X\^2  =  0  (equilibrium)  needed  to  initiate  interface  motion).  Note  that  the  local  driving  force  dG/dt /;  is  concentrated  at  the 
interfaces  because  it  is  zero  for  both  ;/,■  =  0  and  77,  =  1 ;  internal  stresses  are  also  concentrated  at  the  interfaces. 

The  evolution  of  stochastic  initial  data  under  abJ  =  crb2  =  10  GPa  in  the  presence  of  the  equal  square  wave  defect  stress 
fields  ad2  =  adl  =  ±10  GPa  (i.e.,  X^_,2  =  0)  with  b  =  0.2/  and  inclined  at  45°  is  shown  in  Fig.  20.  The  stationary  solution  (last 
row)  consists  of  M,  plates  in  the  regions  of  positive  stress  separated  by  A  plates  in  the  regions  of  negative  stress. 

Instead  of  a  plane  square  wave,  one  can  smooth  the  transition  from  —H  to  H  with  a  continuous  function  defined  in  a  band 
of  width  bt.  In  this  case,  a  change  in  stress  ab2  -  erM  ~  X^_2  below  the  threshold  value  will  induce  only  a  slight  movement  of 
the  interfaces  that  depends  on  the  stress  field  within  bt.  If  b,  <  b,  this  interface  motion  does  not  significantly  change  the 
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is  Mi  and  blue  is  M2.  The  small  islands  within  broad  Mi  plates  are  incompletely  transformed  A  stabilized  by  the  field  Xj^Mi  of  defects.  IfXj^Mi  is  reduced 
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Fig.  15.  Solution  to  the  same  problem  as  in  Fig.  9  but  with  the  periodic  defect  stress  field  =  5cos(167ix)cos(167cy)  and  oA2  =  5 sin (167tx)  sin (167ty)  for 
which  Xj_2  =  -5Aetcos(167r(x  +y)),  a  periodic  function  along  the  diagonal.  Left  column:  ;  second  and  third  columns:  0\  and  cri  -  o2  \  right  column:  the 
local  driving  force,  dG/drj ,,  for  evolution  of  rjx.  The  bottom  row  is  the  stationary  solution  (red  is  Mi  and  blue  is  M2).  (For  interpretation  of  the  references  to 
color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


volume  fraction  of  twins.  But  if  bt  ~  b,  a  variation  in  stress  below  the  threshold  value  changes  the  stationary  microstructure 
and  twin  concentration  continuously  and  significantly.  In  particular,  for  =  5  cos(6nx)  cos(6ny)  and  Cd2  =  5  sin (6nx) 
sin(67iy)  (X^2  =  -5A£tcos(67i(x  +  y)))  (H  =  5  GPa)  and  abi  =  10  GPa,  we  obtained  the  same  stationary  structure  as  for  the 
square  wave  fields.  However,  when  we  varied  <7f,i  between  5  and  15  GPa,  the  concentration  and  widths  of  all  stationary  twins 
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Fig.  16.  Evolution  of  the  initially  stationary  microstructure  obtained  in  the  previous  simulation  (Fig.  14)  after  removal  of  the  defect  stress  field  a&.  At  first, 
for  (Tbi  =  10  GPa,  the  Mi  plates  widen,  but  then,  before  the  transformation  is  complete,  the  boundary  stress  ab i  is  decreased  to  9.9  GPa.  Variant  M2  grows  by 
coalescence  of  the  plates,  finally  converging  to  homogeneous  M2  with  zero  internal  stresses  and  driving  forces.  Left  column:  rjx\  second  and  third  columns: 
o i  and  G\  —  a2\  right  column:  the  local  driving  force  for  evolution  of  r\x,  dG/dr\x.  The  red  corresponds  to  Mi  and  blue  to  M2.  (For  interpretation  of  the 
references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


v  ( mm/s ) 


Fig.  17.  Interface  velocity  v  vs.  the  macroscopic  driving  force  x(  i2  for  four  interfaces  from  the  microstructure  in  Fig.  16. 


varied  continuously  and  coalesced  to  pure  M2  (or  Mi)  at  5  GPa  (or  15  GPa)  (Fig.  21).  We  emphasize  that  the  slow  cyclic 
stress-twin  concentration  curve  does  not  exhibit  a  hysteresis  loop  and  energy  dissipation. 

4.  In  Fig.  22  we  display  the  evolution  of  an  embryo  with  ip  =  0.1  as  in  Fig.  10  but  with  the  addition  of  the  periodic  defect 
fields  <7dl  =  5cos(167tx)cos(167ty)  and  ad2  =  5sm(167tx)sin(167iy).  The  defect  fields  completely  change  the  microstructure 
evolution.  The  stationary  solution  is  a  twinned  structure  with  equal  width  twins,  where  M2  (or  Mi )  is  located  in  regions  with 
positive  (or  negative)  values  ofX,^2  =  -5Aefcos(167i(x  +  y)).  Around  two  opposite  corners  Mi  units  evolve  into  triangles. 

For  the  same  problem  but  with  equal  defect  stresses  erdi  =  ai2  =  10sin{167i(x))sin(167iy)  (X\_2  =  0),  the  symmetry  of  the 
problem  is  changed  from  the  very  beginning,  and  XA  >M  drives  the  microstructure  to  four  twinned  regions,  as  seen  in  Fig.  23. 
One  can  also  see  multiple  nucleation  of  new  M2  (and  Mi)  units  around  the  maxima  of  XdA  >M .  Thus,  the  stress  fields  of  the 
defects  not  only  produce  an  athermal  threshold  but  also  create  sites  for  heterogeneous  nucleation.  However,  the  driving 
force  due  to  decreases  in  the  energy  associated  with  internal  stresses  and  the  surface  energy  exceeds  the  barriers  due  to 
XA  ,M  .  Coalescence  of  Mi  occurs  and  the  final  stationary  solution  is  a  two-twin  microstructure. 

Results  for  the  same  problem  but  with  the  magnitude  of  the  oscillatory  field  increased  to  20  GPa  are  shown  in  Fig.  24.  The 
stationary  microstructure  consists  of  four  twins  with  islands  of  incompletely  transformed  material,  which  may  serve  as 
nucleation  sites  for  austenite  if  it  is  promoted  by  the  boundary  conditions. 
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Fig.  18.  Solution  to  the  same  problem  as  in  Fig.  9  but  with  the  periodic  stress  field  =  od2  =  10sin(87t(x  +y));  hence,  X)  t2  =  0  and  X)  ,A  is  a  periodic 
function  along  the  diagonal.  Left  column: ;  second  and  third  columns:  ct,  and  cr,  -  g2\  right  column:  dC/dr^ ,  the  local  driving  force  for  the  evolution  of  r/, . 
The  bottom  row  is  the  stationary  solution  (red  is  Mi  and  blue  is  M2).  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred 
to  the  web  version  of  this  paper.) 


We  next  considered  a  defect  stress  field  that  is  stepped  in  both  adl  and  ad2  with  ad2  =  -Odi ;  the  defect  field 
=  2Asf(Tc/2  is  plotted  in  Fig.  25.  The  stresses  have  the  values  0  (green),  -10  (yellow),  and  lOGPa  (blue).  The  system 
evolves  into  a  four-grain  microstructure  in  which  many  IVh  variants  from  one  grain  contact  M2  variants  from  other  grains. 

For  the  same  initial  and  boundary  conditions,  we  again  used  a  stepped  defect  field  satisfying  cr^  = -ffdi ,  and 
=  2A£tffd2.  where  the  stresses  assume  the  values  0  (green),  -10  (yellow)  and  10  GPa  (blue),  as  shown  in  Fig.  26;  the 
width  of  each  step  is  b  =  0.2 1.  The  microstructure  evolves  unsymmetrically;  the  stationary  microstructure  consists  primarily 
of  rectangular  lath  and  a  few  needle-like  units  (Fig.  26).  Lath  microstructures  are  usually  observed  when  phase  transforma¬ 
tions  are  accompanied  by  plastic  accommodation  due  to  dislocation  generation  (Ghosh  and  Olson,  1994).  Such  a  structure 
cannot  appear  in  the  absence  of  an  athermal  threshold  because  elastic  energy  minimization  leads  to  relatively  sharp  tips. 
The  widths  of  most  (but  not  ail)  of  the  laths  are  determined  by  the  defect  stress  field. 

Finally,  we  carried  out  a  simulation  similar  to  the  previous  (Fig.  26)  except  that  the  step  width  was  decreased  to  b  =  0.11. 
As  seen  in  Fig.  27,  the  early  stages  of  the  microstructure  evolution  is  similar  to  the  evolution  without  defects  (Fig.  10).  How¬ 
ever,  the  defect  field  arrests  this  complicated  microstructure,  which  is  a  combination  of  four  herring-bone  microstructures. 

In  all  of  the  above  simulations  with  stepped  X^2.  the  initial  stage  of  the  evolution  is  heterogeneous  nucleation  of  M2  (Mi ) 
units  in  the  regions  with  positive  (negative)  Xd  2.  The  late-time  growth  of  martensitic  units  is  driven  by  a  reduction  in  the 
energy  of  the  internal  stress  field. 

Note  that  in  Kartha  et  al.  (1995),  heterogeneous  distributions  of  disorder  and  an  elastic  constant  were  used  to  model  pre- 
martensitic  structures  like  tweeds.  In  Abeyaratne  et  al.  (1996),  an  oscillatory  contribution  was  added  to  a  local  mesoscale 
potential  in  terms  of  the  volume  fraction  of  a  martensitic  variant,  and  this  produces  an  athermal  threshold  for  the  mesoscale 
modeling,  as  in  other  mesoscale  models  (Ghosh  and  Olson,  1994;  Levitas  et  al.,  2004).  This  approach  is  not  applicable  to 
nanoscale  phase  field  modeling  because  it  arrests  unphysical  intermediate  configurations. 


8.  Concluding  remarks 

In  this  paper,  interface  propagation  kinetics  was  incorporated  in  advanced  phase  field  models  developed  in  Levitas 
and  Preston  (2002a, b)  and  Levitas  et  al.  (2003).  For  ID,  we  obtained  and  analyzed  both  analytical  and  numerical 
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Fig.  19.  (top  row)  Stationary  solution  following  evolution  from  a  stochastic  initial  state  under  ob\  =  ab2  =  10  GPa  in  the  presence  of  a  plane  square  wave 
defect  stress  field  (see  text  for  details).  The  final  state  is  a  twinned  microstructure  with  twins  occupying  adjacent  regions  of  positive  and  negative  xj_2. 
(middle  row)  For  ob\  =15  GPa  and  (7b2  =  10  GPa  the  driving  force  for  the  M2  — >  Mi  phase  transformation  exceeds  the  athermal  threshold  of  4  GPa  and 
variant  Mi  grows  until  the  phase  transformation  is  complete,  (bottom  row)  Similarly,  when  ob\  =10  GPa  and  ob2  =  5  GPa,  so  that  the  driving  force  for  the 
Mi  — >  M2  phase  transformation  exceeds  the  athermal  threshold,  variant  M2  grows  until  the  phase  transformation  is  complete. 
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Fig.  20.  Evolution  of  stochastic  initial  data  under  ob\  =  ab2  =  10  GPa  in  the  presence  of  a  plane  square  wave  defect  stress  field  (see  text  for  details).  The 
stationary  solution  (last  row)  consists  of  M,  plates  in  the  regions  of  positive  stress  separated  by  A  plates  in  the  regions  of  negative  stress. 
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Fig.  21.  (a)  Stationary  solution  following  evolution  from  stochastic  initial  data  under  o m  =  ob2  =  lOGPa  in  the  presence  of  the  periodic  stress  field 
den  =  5cos(6nx)cos(6ny)  and  a2d  =  5sin(6nx)sin(6ny )  (H  =  5  GPa);  x}_2  =  -5  Aef  cos  (67r(x  +y)).  The  solution  converges  to  the  same  twinned  microstruc¬ 
ture  as  in  Fig.  19  (top  row),  (b)  Stationary  microstructure  for  crbl  =  14  GPa  and  <Tf,2  =  10GPa.  (c)  Stationary  microstructure  for  abi  =  10  GPa  and 
(jb2  =  6  GPa.  The  stress-twin  concentration  is  continuous  and  the  slow  cyclic  stress-twin  concentration  curve  does  not  exhibit  a  hysteresis  loop  or  energy 
dissipation. 
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Fig.  22.  Evolution  of  an  embryo,  tji  =  0A,  for  <rbl  =  ob2  =  1 5  GPa,  as  in  Fig.  10,  but  with  the  addition  of  the  periodic  defect  field 
( Jen  =  5 cos(16nx) cos(l6ny ),  od2  =  5sinCl6nx)sin(16ny).  The  stationary  solution  is  a  twinned  microstructure  (lowest  row)  with  equal  width  twins,  where 
M2  -  blue  (Mi  -  red)  is  located  in  the  region  with  a  positive  (negative) xj_2.  Around  two  opposite  corners,  Mi  units  evolve  into  triangles.  (For  interpretation 
of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


solutions  of  the  time-dependent  Ginzburg-Landau  equations  for  A-M  and  M-M  interface  propagation.  This  included  the 
determination  of  the  interface  velocity  and  profile,  and  a  study  of  interface  propagation  instabilities  for  thermody¬ 
namic  parameters  near  the  conditions  for  lattice  instability.  In  2D,  the  coupled  Ginzburg-Landau  and  quasi-static 
equations  of  linear  elasticity  were  solved  using  FEM  for  the  evolution  of  microstructures  and  stress  fields  in  square 
samples. 
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Fig.  23.  Evolution  of  an  embryo,  i;,  =  0.1,  for  (7!n  =  ab2  1 5  GPa  (as  in  Figs.  10  and  22)  but  with  the  periodic  field 
=  <Td2  =  1 0 sin(] 6ttx) sin( i 6 Tty)  (Xd  2  =  0).  The  intermediate  microstructure  consists  of  four  twinned  regions  with  clear  boundaries  between  them. 
The  final  stationary  solution  is  a  two-twin  microstructure. 


We  introduced  an  athermal  threshold  in  the  evolution  equation  for  the  order  parameter  and  found  that  it  works  well  in 
ID,  but  this  approach  should  be  checked  in  higher  dimensions. 

We  proposed  the  inclusion  of  oscillatory  background  stress  fields  in  the  phase  field  approach  to  model  athermal 
thresholds  to  interface  motion  due  to  crystal  defects.  With  these  defect  fields  the  system  is  arrested  in  experimentally 
observed  microstructures  instead  of  evolving  to  a  single  phase,  and  rate-independent  stress  and  temperature  hysteresis 
are  present.  In  addition  to  producing  an  athermal  threshold,  these  fields  also  create  sites  for  heterogeneous  nucleation. 
Some  microstructures,  e.g.,  laths,  appear  in  phase  field  simulations  with  oscillatory  defect  fields  that  are  not  seen 
otherwise. 

Similar  results  are  obtained  by  introducing  spatial  oscillations  in  AG"  due,  for  example,  to  fluctuations  in  chemical  com¬ 
position.  An  athermal  threshold  for  twinning  and  dislocation  motion  can  be  introduced  in  the  same  way.  For  other  types  of 
phase  transformations,  oscillatory  electric,  magnetic,  or  other  fields  can  be  used.  The  oscillatory  defect  fields  not  only  stabi¬ 
lize  certain  experimentally  observed  microstructures  and  cause  rate-independent  phase  transformation  hysteresis,  but  they 
also  significantly  change  the  microstructure. 

Oscillatory  defect  fields  not  only  open  up  new  opportunities  for  realistic  microstructure  modeling,  but  also  lead  to  a  new 
major  problem:  how  to  find  a  realistic  oscillatory  field  corresponding  to  a  given  defect  structure.  For  any  given  (for  example, 
experimentally  observed)  distribution  of  dislocations  and  other  defects,  we  can  find  a  numerically  corresponding  stress  field 
<rd  using,  for  example,  the  approach  in  Wang  et  al.  (2001).  For  any  given  (or  evolving)  heterogeneous  alloy  composition,  in 
cases  in  which  the  dependence  of  AG”  on  the  composition  is  known  and  included  in  the  simulation,  the  heterogeneous  (oscil¬ 
latory)  contribution  to  AG"  appears  automatically.  For  a  particular  defect  stress  field,  the  corresponding  fields  of  the  driving 
forces,  Xd^M  and  Xd  _,M .,  point  to  the  possible  effects  of  ad  on  the  microstructure  evolution.  However,  because  the  two  tra¬ 
ditional  contributions  to  the  driving  force  -  the  surface  energy  and  the  energy  of  the  internal  stress  field  due  to  the  heter¬ 
ogeneous  distribution  of  transformation  strain  -  depend  upon  the  evolving  geometry  of  the  martensitic  units,  the  final 
microstructure  may  be  very  sensitive  to  variations  in  Xd^2.  This  is  typical  of  systems  governed  by  material  instability  and 
may  explain  the  variety  of  microstructures  observed  experimentally  that  do  not  follow  from  the  minimization  of  the  Ginz- 
burg-Landau  energy. 

It  is  known  that  the  defect  microstructure  and  phase  transformation  hysteresis  can  be  changed  by  thermomechanical 
treatment  (Hornbogen,  1999);  this  can  be  modeled  with  oscillatory  fields.  In  addition,  if  a  phase  transformation  is 
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Fig.  24.  Evolution  of  an  embryo,  77,  =  0.1,  for  (Tb\  =  ob2  =  1 5  GPa  (as  in  Fig.  23)  but  with  the  amplitude  of  the  periodic  field  doubled: 
(jdl  =  od2  =  20 sin (167tx)  sin (167ty)  (X^2  =  0).  The  stationary  microstructure  consists  of  four  twins  with  islands  of  incompletely  transformed  material, 
which  may  serve  as  nucleation  cites  for  another  twin  if  promoted  by  the  boundary  conditions. 


Fig.  25.  Evolution  of  an  embryo,  7/,  =  0.1,  under  ob\  =  Ow  =  15  GPa  (as  in  Fig.  24)  but  with  a  stepped  defect  stress  field  for  which 
(702  =  —<7id;  X}_2  =  2A£fCrd2  is  shown  top  left.  Stress  values  are  0  (green),  -10  (yellow),  and  10  GPa  (blue).  The  stationary  microstructure  is  shown  in 
the  last  row.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


accompanied  by  plastic  accommodation,  the  defect  structure  changes,  and  the  athermal  threshold  increases  (Ghosh  and  Ol¬ 
son,  1994),  thus  one  needs  to  determine  the  evolving  oscillatory  field.  The  numerical  solution  to  this  coupled  defect  field- 
phase  structure  problem,  is,  however,  quite  challenging. 

In  the  future  we  will  study  oscillatory  defect  fields  in  the  context  of  our  phase  field  model  for  finite  strain  and  material 
rotation  (Levitas  and  Preston,  2005;  Levitas  et  al.,  2009). 
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Fig.  26.  Evolution  of  an  embryo,  r]t  =  0.1,  under  (jftl  =  ob2  =  15  GPa  with  a  stepped  defect  field  satisfying  ad2  =  -0\d  and  Xj_>2  =  2A8t<rd2,  as  in  Fig.  24,  but 
with  a  different  step  distribution  -  compare  X^_2  in  this  figure  and  Fig.  25.  Stress  values  are  0  (green),  -10  (yellow),  and  10  GPa  (blue).  The  last  row  shows 
the  stationary  microstructure,  which  consists  mostly  of  rectangular  laths  and  a  few  needle-like  units.  (For  interpretation  of  the  references  to  color  in  this 
figure  legend,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


X\2  ,  Mi  (T i  CT]-(7"2  Driving  force  for  Mi 


Fig.  27.  Similar  to  Fig.  26  but  with  the  step  width  reduced  to  b  =  0.1  /.  Stress  values  are  0  (green),  -10  (yellow),  and  10  GPa  (blue).  The  last  row  shows  the 
stationary  microstructure,  which  is  a  combination  of  four  herring-bone  microstructures.  (For  interpretation  of  the  references  to  color  in  this  figure  legend, 
the  reader  is  referred  to  the  web  version  of  this  paper.) 
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